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ABSTRACT 

The electric far field radiation pattern is determined for a 
uniformly illuminated, linearly polarized circular aperture 
transmitting through a dielectric radome located in the near 
field of the aperture. A modified electric field integral 
equation is solved using the method of moments procedure and 
the thin shell approximation for dielectrics. The resulting 
solution was computer coded for ogive and spherical radome 
shapes. The program is designed in a modular fashion to 
accommodate the addition of different antenna types, 
illumination functions, or radome shapes. 
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I. INTRODUCTION 


The term radome was coined during World War II as a 
composite of two words, radar and dome. The term originally 
described the class of dome-shaped structures designed to 
house and protect radar ancennas on airborne platforms. 
Current usage of the term has evolved to encompass any 
structure that protects or houses an antenna [Ref. 1]. 

Radomes must satisfy a diverse range of engineering 
requirements. These specifications can be broadly divided into 
two categories: electrical and aero-mechanical. Due to aero- 
mechanical considerations (which are not addressed in this 
thesis) an ogive-shaped radome is commonly found on airborne 
platforms. Chapter II discusses the geometry for an ogive¬ 
shaped body of revolution. Aircraft radomes, especially those 
used at supersonic speeds are subject to mechanical stress and 
aerodynamic heating so severe that the electrical 
characteristics of the radome may be a secondary 
consideration. 

The goal of this thesis is to develop a computer program 
that provides an accurate model for the transmission 
characteristics of a radome, given the geometry and electrical 
characteristics of its material composition. This program is 
called LDBORMM.F, and is described in detail in subsequent 
chapters. 
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A. ELECTRICAL CHARACTERISTICS OF RADOMES 


Electromagnetic transmission characteristics are a primary 
electrical design problem for radomes. In this thesis, the 
radome is modeled as a thin isotropic, homogeneous dielectric 
shell. The reflection coefficient (T) of the radome is given 
by 


r 


_ 3 _ 

2 R g + l\ 


( 1 . 1 ) 


where R s is the surface resistivity of the radome and r] is the 
impedance of free space (120 tt ohms). The surface resistance, 
which can be complex in general, is given by 


R 


e 


£t 

Js 


( 1 . 2 ) 


where E t is the tangential electric field at the surface of 
the radome and J s is the surface current [Ref 2: p. 433]. 

From Equation (1.1), it follows that the reflectivity (T) 
of the radome increases as the surface resistance (R s ) 
decreases. In the limit, where R,=0 (perfect conductor), r=-l 
as expected. When an antenna radia’ -_s through a radome, a 
small portion of the main b r ‘eflected off of the walls 
and forms a lobe. This ' _mmonly called the reflection 
lobe or image lobe. e lobe is usually broader than the 
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The primary concern regarding the image lobe is that large 
discrete reflectors in the detection field can become visible 
to the radar through this side lobe, and be interpreted as 
moving targets by the radar [Ref. 3: pp. 301-311]. The reason 


✓ 

Image lobe / 



In side lobe levels. “V 


Figure 1.1. Radome Reflection Lobe Geometry. 

for false detection of motion from stationary reflectors is 
the apparent doppler of the return of the image lobe since it 
is at a different angle and thus a different frequency than 
the main beam lobe. 


The ideal radome is perfectly transparent to the 
electromagnetic energy that must pass through it. The 






The ideal radome is perfectly transparent to the 
electromagnetic energy that must pass through it. The 
materials available and suited for radome construction in high 
velocity airborne systems are not transparent. Impedance 
mismatches (radome/air interface) and radome geometry alter 
the transmission characteristics of the system. The presence 
of a radome can affect the gain, beamwidth, sidelobe level and 
the direction of the boresight, as well as change the VSWR and 
the antenna noise temperature [Ref. 4: p. 265]. 

B. METHODS OP ANALYSIS 

Most analytic models for radomes are based on microwave 
optics (ray tracing techniques). These models assume a radome 
curvature large enough so that the surface can be considered 
locally flat. In addition, multiple reflections and surface 
waves are usually ignored. These approximations are not valid 
for systems where the radome is in the near field of the 
antenna or the curvature of the radome causes rapid variation 
in the direction of the tangential vector along the surface of 
the radome. Consequently, ray tracing techniques are 
inaccurate for the analysis of high performance antenna and 
radome configurations. Figure 1.2 illustrates the inaccuracies 
inherent to ray tracing techniques. Other analytical methods 
are available, but have various shortcomings and limitations. 
Table 1.1 gives a summary of various analytic methods for the 
solution of electromagnetic scattering problems. 
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The electric field integral equation (EFIE) and the 
method of moments is employed in this thesis. The objective is 
to obtain the solution for the unknown current density which 


TABLE 1.1. SUMMARY OF ANALYTIC METHODS. 


METHOD 

SHAPE " *** 


PPijjpS 

CE- 

Special 

Function 

Series 

Solution 

spheres, 

infinite 

cylinders, 

spheroids 

rigorous 

formulation;separation of 
variable 

techniques ? special 
functions 

Extended 

boundary 

condition 

method 

homogeneous 

objects 

rigorous integral 
equation 

formulation; harmonic 
expansion of Green's 
functions bodies 

Physical 

optics 

large planar 
surfaces 

currents deduced from 
incident fields 

Geometrical 

optics 

large surfaces 

discontinuous field 
patterns obtained at 
shadow and reflection 
boundaries 

Geometrical, 
uniform, 
physical 
theory of 
diffraction 

large convex 
surfaces with 
edges and 
discontinuities 

generalization of 
physical and geometrical 
optics 

Frock theory 

large smooth, 
convex surfaces 

formulation based on 
principle of locality 

Hybrid 

methods 

large class of 
scatterers 

integral equation 
formulation; combines 
method of moments and 
optic derived methods 

Method of 

moments, 

finite 

element and 

finite 

difference 

large class of 
bodies, 
inhomogeneous 
obj ects 

rigorous integral 
equation formulation 








































ray tracing techniques assume 
source rays are parallel and the 
scattering surface is locally flat. 



tangential vectors on the 
surface vary rapidly in 
direction. 


Figure 1.2. Ray Tracing Approximations. 
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occurs in the integral equation. The integral equation is then 
solved using a numerical procedure called the method of 
moments (MM) [Ref. 5: p. 671]. When properly implemented, a 
method of moments solution is considered rigorous; all of the 
scattering mechanisms are included. This thesis adapts a code 
developed by J.R. Mautz and R.F. Harrington for scattering by 
bodies of revolution. 

The proximity of the antenna to the radome requires a near 
field solution of the integro-differential electric field 
equations. Chapter III develops the analytical solution for 
a radome in the near field of a circular antenna and the 
numerical method for its calculation. 

Chapter IV provides details of the program LDBORMM.F. This 
program solves for the far field radiation pattern of a 
circular aperture transmitting through an arbitrary-shaped 
dielectric body of revolution. The EFIE for electromagnetic 
scattering is solved numerically for the system depicted in 
Figure 1.3. 

Chapter V discusses the testing and evaluation of the 
main program. The effects of radome geometry and dielectric 
properties on the far field radiation pattern are examined. 

The appendices provide copies of the source codes and a 
description of variables in the arguments of the subroutines 
and functions contained in the main program. 
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linearly polarized , main beam 
antenna electronically 


scanned 


Figure 1.3 Physical Arrangement for LDBORMM.F. 
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II. OGIVE GEOMETRY 


A. BACKGROUND 

Forward-looking radar systems in missiles and aircraft are 
by necessity located in the nose of the platform. Aerodynamic 
considerations dictate an ogive-shaped radome for most 
applications. The basic problem is illustrated in Figure 1.3. 
In this thesis, the subroutine OGIVE provides the physical 
dimensions of the surface, where the integro-differential 
field equation solutions are satisfied for a loaded (i.e. 
dielectric) body of revolution (BOR). The subroutine 
TESTSPHERE generates a spherical BOR to compare numerical 
solutions of the main program and function subprograms with 
known analytical solutions. TESTSPHERE is used solely for 
validating the subroutines that calculate the circular 
aperture radiation field. It is not used in the calculation of 
the radome scattered fields. In all configurations the antenna 
is centered on the origin of the global coordinate system 
which is the cylindrical coordinate system of the BOR. 

B. DEFINITION OF OGIVE 

The ogive is a figure of revolution obtained by rotating 
an arc of a circle about an axis in the plane of the arc as 
depicted in Figure 2.1. The curve of the ogive is rotated 
around the z axis to produce a surface of revolution. 
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The profile of the surface is generated in the r-z plane 


by the equation 


r(z) •jR 2 -z i *b-R 


( 2 . 1 ) 



OGIVE body of revolution (bor), r(z) is rotated 
around z axis The radius (R) of the parent 
circle originates in the z ■ 0 plane. 


Figure 2.1 Ogive as Defined by Equation (2.1). 


where R is the radius of the parent circle and b is the radius 


of the base. 





C. MODIFIED OGIVE FOR LDBORMM.F 

In the program LDBORMM.F the antenna is located at the 
origin of the coordinate system. In a real system, however, 
the antenna may be displaced on the z axis, from the origin of 
the defining Equation (2.1). The subroutine OGIVE prompts the 
user to enter the value z', which is the displacement on the 
z axis of the antenna location and the position on the z axis 
where the radius (R) of the parent circle of the ogive 
commences rotation. Figure 2.2 depicts the modified profile of 
an ogive as calculated by the subroutine OGIVE. 



Figure 2.2 BOR Generated by Equation (2.2). 


The modified definition for the radius as a function of 
position on the z axis is given by, 
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r{z) 


b 


( 2 . 2 ) 


zzz' 

Jr 2 -(z-z') 2 +b-R z>z' 

where z' is the displacement on the z axis. Using this 
modified definition the source was located at the origin of 
the cylindrical coordinate system used to generate ogive(i.e., 
the source is located in plane z=0). To reiterate, the z' of 
Equation (2.2) and Figure 2.2 denotes the distance from the 
origin where the radius of the parent circle is defined and is 
not related to the source coordinates. 

D. SUBROUTINE OGIVE 

The subroutine OGIVE returns values for the global 
variables NP, ZH and RH. NP is the number of surface 
generating points on the ogive in the r-z plane. In the 
subroutine the range of NP is limited to 3 < NP < 400. The 
program prompts the user for the following inputs: radius of 
the parent circle (R), z', and the base radius (b). 

A standard rule of thumb for the convergence of the series 
approximation for the surface current is that the spacing of 
the surface generating points be approximately one-tenth 
wavelength or less. The following algorithm computes the 
number of points required for a segment size of approximately 
one-tenth of a wavelength on the BOR 
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Z = \J2 *R*b-b z 
0 = ARCSIN 1 

(2.3) 

AL = B*R + Z' 

NP = 2 * INTEGER (5 *AL) +1 . 

NP is the global variable for the number of points in the r-z 
plane of the ogive [reference Figure 2.2]. For a circular 
aperture, the field incident on the radome tends to vary more 
rapidly as the angle off boresight increases, thus more 
generating points are needed. By slightly modifying the above 
algorithm, a surface with any arbitrary segment size can be 
generated. 

Equation (2.3) gives the algorithm employed by the 
subroutine OGIVE to calculate the number of points (NP) on the 
ogive surface in the r-z plane. Data generated uniformly along 
the radome surface will be nonuniform in z, r, and 0 . The 
subroutine TESTSPHERE generates data points with a uniform 
spacing in angle 0 , whereas OGIVE points are uniform along S. 
Thus a comparison of the two outputs will not yield coincident 
data points. Appendix A contains the source code for 
subroutines OGIVE and TESTSPHERE. Appendix E provides 
mathematical details relevant to the distribution of data 
points in angle 0 for the subroutines OGIVE and TESTSPHERE. 
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If the number of points is greater than 400, array limits 
are exceeded and the main program will write a warning to the 
screen and stop the program. For bodies of revolution where 
the arc length is greater than 40 \ the dimensions of the 
following global variables must be changed accordingly 
throughout the program: RH, ZH, Z, R, B, C, ZLO and ZL. ZH n 
and RH n are the product of the local coordinates z n , r n 
(computed by the BOR subroutine) and the wave number k=27r/A.. 
All spatial dimensions throughout the main program, 
subroutines and functions are in wavelengths (A). 

For surface shapes with definitive geometries other than 
an ogive or sphere the user must provide an appropriate 
subroutine. The subroutine must return the global variables 
NP, ZH, RH, b, R and Z'. These variables are required to 
define the generating points of the BOR profile (see Eqn.2.2). 
The call for the subroutine must be inserted in the logical 
block of code in the main program which establishes the BOR 
geometry. The subroutine is appended to the main program. 
Chapter V addresses modification of the main program for 
various bodies of revolution. The program LDBORMM.F saves the 
inputs from the BOR subroutine and writes them to the 
external sequential file OUTLDBOR. 
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III. METHOD OF MOMENTS 


A. OVERVIEW 

In order to calculate the scattered field pattern of a 
dielectric radome the electric current (J s ) must be integrated 
over the surface of the body of revolution. The magnetic 
vector potential (A) and scalar electric potential ($) are 
defined in terms of the surface current 

££(*,0,$) = -jo>4 (Jg) - (3.1) 

where E s (R,0,0) is the scattered electric field in spherical 
coordinates, and L(») is an operator introduced for notational 
convenience. The potentials are 

dS> ( 3 * 2 ) 

s 


0 ITS ds ' ( 3 . 3 ) 

s 

0 = —J- V e . (3.4) 

R is the magnitude of the difference r and £',the position 
vectors of the field and source points respectively. The 
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operator ?,• is the surface divergence on S [Ref. 6: pp. 15 - 
16]. Figure 3.1 illustrates the quantities. 


P(r,e,$> 



To calculate the electric Held , Cue to a source 
with an electric current density J, 

Equations (3. 1 ) through (3 4) must be calculated. 

Figure 3.1. Scattering due to Surface Currents. 

The objective is to solve for the current on the surface 
of the radome. Following the method of moments (MM), the 
electric current (J t ) on the surface (S) is represented by 

iZs = 5Taj(^ J JJ iZol + * Zai) • (3.5) 

in Equation (3.5) and are the tangential and 

azimuthal components of the current on the surface of the 
radome. The expansion functions for are chosen to be 


L6 








pulse function 

Figure 3.2 Basis Functions Pulse (P) and Triangle (T) of 
Equations (3.6) and (3.7) 
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j‘=£ I£B. e Jn * 7 = 1,2,- NP-2 (3.6) 

r i 

= $ JW j = 1,2- - NP-1 (3.7) 

r j 

where the azimuthal mode number is n=0,±l,±2. .. .±«>. An exact 
solution requires an infinite number of terms in the sum over 
the index n. In practice, when the source (antenna) is on the 
axis of symmetry, the sum converges rapidly and only a few 
azimuthal modes are required. 

The unit vectors t and $ are in the tangential and 
azimuthal directions. The functions Tj(t) and Pj(t) are the 
triangle and pulse functions as shown in Figure 3.2. The 
abscissa t is the arc length along the generating curve of the 
BOR. It is assumed that the generating curve consists of NP-1 
straight line segments where NP is an integer. The expansion 
functions of Equations (3.6) and (3.7) are especially 
appropriate if the BOR is an infinitely thin perfectly 

A 

conducting surface with edges. This is true because the t 
directed electric current approaches zero at an edge, whereas 
the 0 directed electric current may grow large [Ref. 7: p. 6]. 
The MM solution of the EFIE reduces the integral equation to 
a matrix equation for the unknown coefficients I 1, and I* . This 
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is accomplished by the MM testing procedure. The test 
functions (W^,^/) are chosen to be the complex conjugates of 
Equations (3.6) and (3.7). Equation (3.5) is substituted 
into Equations (3.2) to (3.4) for the solution of Equation 
(3.1). The dot product of Equation (3.1) is then taken with 
the test functions These dot products are then 
integrated over the surface S. The resulting matrix equation 
is 


■j 



z *♦' 

-1 

V e ' 

J* 


Z*t 

z 




where Z is a square matrix called the impedance matrix, and V 
is a column vector called the excitation vector [Ref. 7: pp. 
6 - 31 ] . 

The impedance elements are identical to those described by 
Mautz and Harrington. Detailed equations for them appear in 
reference [Ref. 7] and will not be repeated here. However, the 
excitation elements will depend on the characteristics of the 
antenna and the location of the radome. The derivation of the 
excitation elements for a BOR in the near field of a source 
has not been presented elsewhere, and therefore will be 
derived in detail in section (3.B.). 

In LDBORMM.F, the subroutine ZMAT calculates the impedance 
matrix and the subroutine GENEX calculates the excitation 
vector (V). The subroutines DECOMP and SOLVE solve the system 
of equations resulting in the electric current coefficients 
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(It,*) . T ^ e curre nt coefficients are substituted in Equation 
(3.2), which determines the surface current (J s ) [Ref. 7: 
pp.41-64] . 

With the surface current (J.) calculated, the scattered 
electric field (E s ) is determined by Equation (3.1). The total 
electric field, E(R,0,0), is then given by 

£(J?,6,4>) = El * El < 3 - 9 > 

where E a is the electric field of the circular aperture. 

B. NEAR FIELD OF A CIRCULAR APERTURE 

In order to maintain a clear picture of the analysis that 
follows, it is advantageous to keep in mind the following: 

1. The electric current distribution over the circular 
aperture, which represents an ante^? _, is defined in Cartesian 
coordinates. 

2. The electric field due to the aperture and the scattered 
electric field due to the radome are computed in spherical 
coordinates (R,0,<£). 

3. Since the ogive is a body of revolution it is generated in 
a cylindrical coordinate system (r,0,z). 

4. Primed quantities denote source coordinates, unprimed 
quantities denote coordinates at the point of observation 
(with the exception of z', which was discussed earlier). 

To calculate the excitation vector (V) of Equation (3.8), 
the electric field of the antenna must be known. Since the 
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antenna is typically operating in the near field of the 
radome, far field approximations are not valid and a more 
accurate method of calculation must be employed. Since many of 
the radar antennas used in airborne applications are circular 
paraboloids or slotted waveguide plates, the radiation source 
is modeled as a circular aperture with a uniform current 
distribution polarized in the x direction 

Jix'.y'.z') = J x X . (3.10) 

For simplicity the electric current distribution over the 
aperture is assumed to have constant amplitude and a linear 
phase 

j x = j Q e jkx ' B ine.coB*. ( (3.11) 

J 0 is the amplitude of the electric current, while 0 S and 
are scan angles in the spherical coordinate system. The 
coordinates of a point in the plane of the aperture are given 
by (x',y',0). 

To compute the Cartesian components of the electric field, 
J x is used in Equations (6-108a) to (6-108c) of Reference 
[Ref. 5: p. 283]. A straightforward transformation of 
rectangular to polar coordinates (see Appendix E) yields the 
following expressions for the Cartesian components of the 
scattered electric field 
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where 


E = ~ Jr] J ° 
* 4nJc 


JJe XJ [G 1 + (r cos^-r' cos^G^i' dz' d$’ (3.12) 


F _ -JiK 

y 4JtJc 


^^2 j J e xl (z cosfy-z' costy) 

(z sin^-r' sin^O^r 7 dz'dfy' 


(3.13) 


_ -j'n^o 


jJ e xl (z cosfy-z' cosfy) z G 2 z' dz' d$' (3.14) 


xl = jk (r' cos<J) / sin0 fl -i?) 


(3.15) 


R = J (r -r') 2 + z 2 +4 rr 




(3.16) 


Q1 = k 2 R 2 -l-jkR 
R 3 


(3.17) 


G2 = 


3 +3jkR~k 2 R 2 
R 5 


(3.18) 


Figure 3.3 defines the quantities used in the equations. It is 
assumed that the antenna is scanned such that <p s = 0. 

The rectangular components of the electric field are 
converted to spherical components using the transformation 
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Eq ■ E^cosOcos^+ £ y cos0sin$-E f sin0 


( 3 . 19 ) 


E k - £ x sin0cos$-*-E y sin0sin$ + E t cos0 


- -£ x sin$ + E y C 08 $ . 



The function subprograms 

CIRCTHETA , CIRCRHO and CiRCPHi calculate the spherical 
electric field components at the point ( P ) 

The point (P) Is defined in cylindrical coordinates, 
r(z). 0 and Z 

The primed coordinates are In the plane Z • 0 
The circular aperture Is centered at the origin 

Figure 3.3 Geometry for Equations (3.12) to (3.18) 

The function subprograms CIRCTHETA, CIRCRHO, and 
calculate the spherical electric field components E e , 
E # , respectively. 


(3.20) 

(3.21) 


CIRCPHI 
E r , and 








1. FUNCTION SUBPROGRAMS, CIRCTHBTA, CIRCRHO, AND CIRCPHI 

The function subprograms CIRCTHETA, CIRCRHO and 
CIRCPHI are called in the subroutine GENEX. In the main 
program, the vector field due to the source (E a ) is added to 
the scattered field vector (E“) of the body of revolution to 
obtain the total field (E R e *) , per Equation (3.9). Since the 
far field approximation is valid in this case, the closed form 
solution for a circular aperture derived in Section III.B.l.b 
can be used. Subroutine GENEX requires the components of E a in 
the near field. The function subprograms calculate the 
spherical electric field components due to a uniformly 
illuminated circular aperture at an arbitrary point in space 
specified by the spherical coordinates R, 0, 0. The source can 
be scanned off boresight by entering non-zero scan angles when 
the program is executed. 

In the subroutine GENEX, these functions are called 
for each point of integration in the tangential and azimuthal 
components of the excitation vector (V) which occur in 
Equation (3.8), and are described in the following section. 

Appendix D contains a table of the variables in the 
argument list of the function subprograms and the source code 
for CIRCSUB.F. The program CIRCSUB.F was used to validate 
the algorithms of CIRCTHETA, CIRCRHO and CIRCPHI. The 
programs TESTCIRC.M and SCANPLT.M (Appendix D) were used to 
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generate the analytic solution and plot the comparative 
results. 

a. Numerical Method 

To determine the electric field of the antenna and 
the excitation vector requires evaluating several integrals. 
In general these integrals can not be reduced to a closed form 
expression for the geometries under consideration. The 
integrals must therefore be integrated numerically. The chosen 
method for evaluating all integrals in this program is 
Gaussian quadrature. In Gaussian quadrature algorithms, the 
integral is over a normalized interval, such that -1< x < 1. 
N-point Gaussian quadrature approximates an integral by 

/„ * EL ». 13 • 221 

where 


N = 2/n + 2 (3.23) 

and are the coefficients and Xn are the zeros of the 
Legendre polynomial of order m (m=0,+l,±2.. .<») [Ref. 8]. 

The number of terms in the sum is always an even 
number. The program GAUS.F (see Appendix C.) will generate an 
external sequential file for any range of N where 0 < N < 500. 
The first record in the file is an even number (N) . The 
remainder of the file contains N rows of x n and H • The gaus 
files are stored in a subdirectory named "gaus" of the 
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directory containing the main program LDBORMM.F. The path to 
the data files is "(MAIN DIRECTORY)/gaus/gaus###, where ### is 
N. To create the data file (gaus###), execute GAUS.F in the 
subdirectory "gaus." The screen prompts the user for the data 
filename and an even number of points. 

For integrals that are over a range other than -1 
to +1, a change of variables is necessary before Gaussian 
quadrature is applied. Integrations over the antenna aperture 
(from zero radius to the antenna radius r c ) require this 
change of variables is given by [Ref 8: pp. 523-527] 



A similar technique is used to evaluate the <t> 
integrals. The interval of integration for 0 is symmetric and 
does not require a change of variables since, -jt< <p < n 

. (3.25) 

Equations (3.12) through (3.14) are integrated by this method. 
In the function subprograms, CIRCTHETA, CIRCRHO and CIRCPHI, 
the variables of integration over the aperture are PHIPRIME 
( 0 ') and RP (r') and the limits are from 0 to 2 ir and 0 to r 0 , 
respectively. 
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In order to verify the algorithms and coding for 
CIRCTHETA, CIRCRHO and CIRCPHI the field magnitude E e was 
computed and compared with the analytical solution given by 


jkr a 2 E 0 e~ jk * J x (icr 0 (sin0-sin0 5 ) 
r kr Q {sinO- sin0 # ) 


COS (0) (3.26) 


where E 0 =r?J o is the electric field amplitude, is the first 
order Bessel function, 0, is the scan angle and (R) is the 
distance from the center of the aperture to the far field 
point (P) of Figure 3.3. It is assumed that R is much greater 
than r' in the above formula. Therefore any vector from a 
source point to the observation point is approximately 
parallel to the vector R [Ref. 9: pp. 478-487]. 

The far field distance for which Equation (3.26) is 
valid depends on radius (r 0 ) of the aperture. Antenna 
engineers commonly consider the far field to exist at 
distances greater than 



However in optics the far zone applies to values of R where 
the Fresnel number is less than one. Using the latter 
criterion the minimum distance for which Equation (3.26) is 
valid is given by [Ref. 10: pp. 95-96] 
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(3.28) 


1 i 


kr 0 2 \ 2 


I”. 


2nR 

2 


The program TESTCIRC.M in Appendix D is used to 
calculate the analytical solution for E 9 in accordance with 
Equation (3.26). The results are compared with the calculated 
field computed by CIRCSUB.F. In Figures (3.4), (3.6), and 
(3.8) the analytical solution is plotted with a solid line 
while the values calculated by the test program CIRCSUB.F 
are plotted as points. 

Figures (3.4), (3.6), and (3.8) show the 
convergence of the numerical solution from CIRCSUB.F as a 
function of the number of integration points with the 
analytical solution given by Equation (3.26) for one, two and 
five wavelength radius circular apertures. Figures (3.5), 
(3.7), and (3.9) show plots of the field magnitudes (E 9 , E R , 
E t ) versus the angle in degrees. Figures 3.10 through 3.12 
show plots from CIRCTHETA and Equation (3.26) for various scan 
angles and aperture radii. (Note that the vertical scale in 
Figures 3.4 to 3.12 is in volts/meter, not dB.) 

Table 3.1 summarizes parameters for Figures (3.4) 
through (3.12). The plotted results verify the excellent 
agreement between the algorithm of the function subprograms 
and the analytical solution. It is important to note the 
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number of points of integration in r' (CNRHO) and <p' (CNPHI) 
of Equations (3.13) through Equation (3.15) required for 
convergence to the analytical solution of Equation (3.26). 
The number of points of integration in the table are the 
minimum number required for convergence. As the radius of the 
aperture increases, the number of iterations required 
increases as the product of CNRHO and CNPHI. For example, 
Test 3 requires 1200 iterations for convergence to the 
analytic solution. In the program LDBORMM.F the functions 
TABLE 3.1. DATA SUMMARY FOR TEST CASES. 



subprograms are called five times in the subroutine GENEX 
which is in turn called once for each azimuthal mode. Thus the 


number of calls for the function subprograms is the product of 
the number of integration points in t (NT), number of 
integration points in <p (NPHI) , the number of azimuthal modes 
(MODES) and the number of generating points on the surface of 
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the BOR. For a large BOR, the computer run time can become 
significant. It is therefore prudent to test the function 
subprograms for the minimum number of points of integration 
required for convergence before executing the main program. 
2. THE SUBROUTINE GENEX 

The subroutine GENEX calculates the excitation vector 
(V) for Equation (3.8) . All variables in the argument list are 
inputs with the exception of the excitation vector. Appendix 
B contains a list of variables in the argument of the 
subroutine GENEX. 

The jth element of V 1 and V* are given by 

Vn P = - (f Win 'El ds, j-1,2, . . .NP-2 (3.29) 

^' p = ![ s ^n' El ds, j-1,2, . . . NP~ 1 (3.30) 

where the indices j and n are the same as Equations (3.6) and 

(3.7) , and W L * is the testing function. The test functions are 
chosen by Galerkin's method and therefore are the complex 
conjugates of the expansion functions Q} 4 (Equations (3.6) and 

(3.7) ]. E p a is the pth component electric field of the 
antenna and S denotes the surface of the radome. The 
representation p is a generalized for the spherical components 
R, 0 and <p . These components are provided by the subroutines 
CIRCTHETA, CIRCRHO, and CIRCPHI. 
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Figure 3.4 Test 1. 
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Figure 3.6 Test 2. 
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Figure 3.7 Test 2, E(R,6,$) Calculated by CIRCSUB.F. 
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Figure 3.9 Test 3, E(R,0,0) calculated by CIRCSUB.F. 
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TEST OF SUBROUTINE CIRCTHETA 



Figure 3.10 Test 4, 0 Scanned to 30 Degrees. 
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TEST OF SUBROUTINE CIRCTHETA 
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Figure 3.11 Test 5, 0 Scanned to 45 Degrees. 
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TEST OF SUBROUTINE CIRCTHETA 



Figure 3.12 Test 6, 6 Scanned to 60 Degrees. 
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The quantities E R/ E s and E # completely define £‘ and are used 
by GENEX to evaluate Equations (3.28) and (3.29) numerically. 
To further reduce equations (3.28) and (3.29) to a form 
suitable for programming requires evaluation of the dot 
products. Referring to Figure 3.13, the tangent direction 
along the surface can be written as 

£ = f sinVj + 2 cosvj (3.31) 

By the orthogonality properties 

$ • £ * 0 

| •0 = 0 (3.32) 

t • $ = 0 . 

Therefore, from (3.28) and (3.29) v t¥ =v w =v^=0. The remaining 
elements are 


Vai‘ f e' jn *\ f Ti CF x E r dt + f TJ CF U1 E r at Id* (3.33) 

[■’A, •'*1.1 

Vni = f e f T* SF s Eq dt+[ T] SF^ Et dtl d<|) (3.34) 

v*f=[ e- jn * f 5 P i dt d$ (3.35) 

JAj ? J 

A, = t;.i - t; (3.36) 
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^1*1 = ti* 2 ~ 1 


( 3 . 37 ) 



SF l = sin (- 00 (3.39) 

CF X = COS (V i - 0') • (3.40) 

The quantities are defined in Figures 3.2 and 3.13. The 



Figure 3.13 Geometry for Equations (3.30) through 
(3.39) 


elements of the excitation vector (V) are stored in the global 
variable R, where 
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V i t9 +V i tR is in R(i) i=l to NP-2 

V ** is in R(i+MT) i=l to NP-l 

and MT=NP-2 [Ref. 7: p. 57]. 

C. MODIFICATION OF THE IMPEDANCE MATRIX FOR DIELECTRICS 
1. THIN-SHELL APPROXIMATIONS FOR DIELECTRICS 

For a perfect conductor the surface resistance (R s ) is zero. 
To satisfy the boundary condition requires 

^ = -E* (3.41) 

for the tangential components of the incident and scattered 
fields. In this case the incident field is that due to the 
antenna. When a dielectric body is present, the incident 
electric field produces a polarization current rather than a 
conduction current. In this case the boundary conditions state 
that the electric and magnetic fields must be continuous. The 
total electric and magnetic fields are 


E = El + El 
H - HI * HI ■ 

Maxwell's equations must be satisfied where 

V x £ = -ju\i 0 H 


(3.42) 


(3.43) 


and by superposition 
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v X K1 = -ju\L 0 HI 


( 3 . 44 ) 


so that 


V x H = jo>e 0 E + jo) (e -e 0 ) E (3.45) 

where the second term on the right side of Equation (3.44) is 
nonzero for e*e Q . 

The polarization current may be defined as 

i7p * jo) (e ~t 0 ) £ . (3.46) 

Maxwells' second equation becomes 

V x H* = + Jp . (3.47) 


Thus E s and H s can be determined from the polarization current 
(J p ) radiating in free space. For dielectric bodies with e r >>l 
or very thin dielectric shells with thickness t << A, the 
tangential component of the polarization current will be much 
greater than the normal component. In these cases the normal 
component can be neglected and a modified form of the EFIE 
results 


L(J P ) +■ 


El 


(3.48) 


JO) (e -e 0 ) t - 

where the coefficient of J p is called the surface loading or 
surface impedance (R s ) . (Frequently it is referred to as a 
surface impedance because it can be complex. Here it will be 
called a resistance to avoid confusion with the surface 
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impedance approximation in which both electric and magnetic 
currents are allowed.) For a lossless material R s is 
imaginary; for a resistive material it is real. 

The new matrix equation corresponding to Equation 
(3.8) is 


[J] = [Z m + Z h )-' [V\ (3.49) 

where 


Z m = Jf s ita'UJj ds (3.50) 

= ffjL ■ RsJ^ds . (3.51) 


2. EVALUATING THE IMPEDANCE MATRIX 

The impedance matrix (Z) of Equation (3.8) is 
calculated by the subroutines ZMAT, ZLOAD, and ZTOT. If for 
some reason the impedance of the radome is zero (a perfect 
conductor), only the original impedance matrix from ZMAT is 
needed. For other materials the load impedance is calculated 
by ZLOAD and added to the ZMAT result the subroutine ZTOT. 

The thin shell approximation assumes the thickness(t) 
is much less than the wavelength (A) and the skin depth (6) of 
the shell material. This is a valid assumption for most 
practical radome designs and materials. The polarization 
current is primarily tangential to the surface of the 
dielectric and the normal components are negligible. This 
approximation is applicable to lossy as well as loss free 
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dielectrics. In the lossy case the load matrix (Z L ) has a real 
and imaginary component, while the loss free case is purely 
reactive. 

At this point it is necessary to relate the surface 
impedance to the electrical properties of the material such as 
the dielectric constant and loss tangent. The electric field 
due to a thin dielectric shell is given by equation (3.48) 
with R s defined by 


r = _—±- 

jto A e t 


(3.52) 


where Ae=e-e 0 [Ref. 11: pp. 531,532]. As stated earlier, the 
dielectric constant of the shell may be complex. Letting 
e=e'-je" gives 


*5 = 


where e ' =e n e. 


ju (t'-jt"-t 0 ) t 
The electric loss tangent is given by 


(3.53) 


tan (6) = = —— (3.54) 
t 

where <5 is the loss angle and a=ue H is the conductivity 
representing all losses in the medium. Therefore the loss 
tangent is a measure of the power loss in the medium. A medium 
is a good conductor if a» coe, and a good insulator if a<<coe 
[Ref. 2: pp. 342,343]. For <oe 0 * 1/60A, and thickness t=nA 
substitution of (3.58) into (3.53) yields 
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(3.55) 


s n [ j (e r -l) +e r tanfi] 

Typical radome materials have 2 < e r < 10 and 0.0005 < tan <5 
< 0.005 . Inhomogeneous radomes (i.e. sandwiches) can be 
treated by using an equivalent dielectric constant. 
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IV. PROGRAM DESCRIPTION 


A. DATA FLOW 

The program LDBORMM.F has eight subroutines and four 
function subprograms. Subroutines and functions not covered 
here are documented in reference [Ref. 7]. Table 4.1 gives a 
brief description of all subroutines and function subprograms 
in LDBORMM.F. 


TABLE 4.1 SUMMARY OF SUBPROGRAMS 


NAME 

TYPE 

DESCRIPTION 

OGIVE 

subroutine 

defines the radome dimensions 

GENEX 

subroutine 

calculates the excitation vector 

ZMAT 

subroutine 

calculates the moment matrix 

ZLOAD 

subroutine 

calculates the submatrix for the 
dielectric BOR 

ZTOT 

subroutine 

calculates the moment matrix for 
dielectric BOR 

DECOMP 

subroutine 

inverts the moment matrix 

SOLVE 

subroutine 

calculates the current 
coefficients 

PLANE 

subroutine 

calculate E far field 

BLOG 

function 

fourth order series expansion of 
log(x), for x< 0.1. 

CIRCTHETA 

CIRCRHO 

CIRCPHI 

functions 

calculate the p, 0 and <p 
components of the electric field 
at a point (P). 







































The block diagram of Figure 4.1 shows the structure of the 


calls for subroutines and functions in the main program. 



Figure 4.1 Structure of main program LDBORMM.F. 
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B. PROGRAM PARAMETERS 


1. NUMBER OF POINTS ON BOR 

The number of points (NP) must be at least 3 but less 
than or equal to 400. This restriction is due to the array 
dimensions set in the code and is not a limitation in the 
solution. The subroutine OGIVE generates segments on the 
surface of the BOR at approximately one tenth of one 
wavelength. Thus, a BOR with an arc length greater than 40 X, 
matrix dimensions are exceeded for the global variables RH, 
ZH, Z, R, B, C, ZLO, and ZL. The subroutine TESTSPHERE allows 
the user to select the number of points on the surface of a 
sphere. The points generated by TESTSPHERE are spaced 
uniformly in the angle 0. 

2. DIRECTORY STRUCTURE 

a. Subdirectory GAUS 

A subdirectory named "gaus" must contai i the 
external sequential files with the weights and abscissas for 
the Gaussian quadrature algorithm. The data files are 
generated by the program GAUS.F, which is executed in the 
subdirectory "gaus." Appendix C contains the program GAUS.F. 
The status of the files in the main program is "OLD", thereby 
invoking an error condition if the file does not exist. Figure 
4.2 illustrates the directory structure. 
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External Input/Output 


The external sequential data files required for the 
main program are shown in Figure 4.2. The programs that 
generate the files are indicated by dashed lines and arrows. 
The dashed line originates at the generating program and the 
arrow terminates on the file generated. 

c. The Program TEST.M 

The programs TEST.M and TESTCIRC.M must be executed 
before the plot routines GENPLT.M and SCANPLT.M. The program 
TEST.M (or TESTCIRC.M) generates the analytical solution given 
by Equation (3.26). The parameters for BOR radius, antenna 
radius and scan angle are entered by the user . These 
parameters must be the same as the inputs to the main program. 
The routine GENPLT.M plots the analytical solution generated 
by TEST.M and the numerical solution calculated by the main 
program for comparison. A comparison of the two results 
illustrates the effect of the radome on the antenna pattern. 

d. The Subdirectory CIRCTEST 

The subdirectory CIRCTEST contains the programs and 
data files required to test for the number of terms necessary 
for the series approximation of the EFIE to converge to the 
analytic solution. The main program requires extensive 
execution time for large BOR or large antenna. In order to 
minimize run time, the minimum number of integration points 
over the antenna is determined by the execution of the 





programs in the subdirectory CIRCTEST. The plotting routine 
SCANPLT.M plots the results for comparison. 

Table 4.2 provides a summary of the files in this thesis. 


TABLE 4.2 FILE 

SUMMARY 




ldbormm.f 

main dir./ 

main program, solves for 
the radiation pattern of 
radome and antenna 

outldbor 

main dir./ 

external data file for 
output Of LDBORMM.F 

bcpole.m 
bxpole.m 
bang.m 

main dir./ 

data generated by 

LDBORMM.F for plotting 
with GENPLT.M 

test.m 

main dir./ 

program to generate 
solution for Eqn.(3.26) 

ctest.mat 

main dir./ 

data generated by TEST.M 

genplt.m 

main dir./ 

plotting routine 

gaus. f 

main dir./ 
gaus/ 

generates gaus### files 

gaus### 

main dir./ 
gaus/ 

. 

data files for Gaussian 
quadrature integration 

circsub.f 

main dir./ 
circtest/ 

program to test function 
subprograms for minimum 
number of integration 
points 

testcirc.m 

main dir./ 
circtest 

generates - ration for 

E C"■ ' ‘ . 

circ.mat 

main dir./ .rated by 

circtest ,.:.M 

scanplt.m 

main dir./ 
circtest 

plotting routine 


C. MAIM PROGRAM MODIFICATIONS 

The main program is designed in a modular fashion in order 
to facilitate modifications for other radome shapes, antenna 
types or dielectric profiles. 
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1. Modifications for Radome Shapes 

The main program contains subprograms to generate 
ogive-shaped or spherical bodies of revolution. When executed 
the main program prompts the user for the type of BOR desired. 
Additional types of BORs may easily be added to the existing 
menu. The following modifications are required: 

(a) Insert code in main program to write selection number and 
type of shape to screen menu. 

(b) Insert code in logical block that calls subroutine 
generating selected shape. 

(c) Append subroutine program to main program. Subroutine must 
return the variables NP, RH, ZH, b, R and Z'. 

These variables are required for subsequent calculations and 
data entry to external files. 

2. Modifications for Antenna Types 

To modify the main program for an antenna type that is 
not a circular aperture: 

(a) Substitute the alternate function name in the main program 
where ETF and EPF are assigned. These are the far field 0 and 
0 components of the antenna, so a closed form approximation 
can be used if one exists. 

(b) Substitute the alternate functions in the subroutine 
GENEX, to calculate the variables SI, S2, S3, S4 and S5. 

(c) Append the functions to the main program. 
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To be consistent with the model developed earlier, the 
function subprograms must calculate £(R,0,0) in spherical 
components at points in the near field. 

3. Modification of Dielectric Profile 

The vector ZLO has a number of elements equal to the 
number of segments (NP-1) on the surface of the loaded body of 
revolution. The elements in the array are the complex values 
of the surface resistance (RJ on each segment of the BOR. As 
written, the main program prompts the user to enter the 
complex surface impedance when the main program is executed. 
As presently coded, the entered value is assigned to each 
segment. In order to examine the effects of a nonuniform 
surface impedance: 

(a) Modify the program to read an external sequential file of 
length (NP-1), with the values of the desired impedance 
profile. 

(b) Store the values in the array ZLO in the main program 
before calls to ZLOAD and ZTOT. The order of storage in ZLO 
must correspond to the order of the coordinates stored in RH 
and ZH. 
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V. SUMMARY OF RE8ULTS 


A. VALIDATION OF PROGRAM 

In order to validate and test the main program, the 
surface resistance was assigned a very large value. For very 
large values of surface resistance, the scattered electric 
field (E s ) approaches zero. The radiation pattern for a BOR in 
the far field with very large surface impedance is therefore 
the far field pattern of the antenna. For a circular aperture, 
the far field radiation pattern is given by Equation (3.26). 

Figures 5.1 through 5.3 compare the values for the 
radiation pattern calculated by the main program and the 
analytic solution generated by TEST.M. The calculated values 
are plotted as (.) and the closed form solution to (3.26) is 
a solid line. The ordinate axis units are volts/meter. The 
abscissa is the angle from the axis of symmetry 0 in degrees. 
The plots of Figures 5.1 to 5.3 show the magnitude of the 
electric field. The actual field values are complex and it is 
the complex values that are used to calculate the excitation 
elements. The figures show the parameters NP, NT, NPHI, CNRHO, 
CNPHI, SCAN, and ARAD as defined in Appendix B. The angle <p 
of the receiver is zero for all test cases. The number of 


55 



















F LDBORMM 



58 












points of integration for the antenna were determined by the 
program JIRCSUB.F to insure that aperture integrations 
converged. Figures 3.4, 3.6, and 3.8 show the results. The 
number of points of integration needed over the BOR surface 
(NT and NPHI) becomes very large as the radial extent (hence 
circumference) of the radome increases. In addition to an 
increase in the number of integration points required, the 
number of azimuthal modes also increases for a large BOR. In 
all test cases the mode number is one (n=±l). Test results 
indicate more points of integration and higher modes are 
required for the calculated values to converge to the analytic 
solution as the antenna radius or sphere radius is increased. 
The execution times were in the range of 2 to 4 hours on a Sun 
Sparcstation 2. 

Figures 5.4 through 5.6 show plots of the magnitude of the 
spherical electric field components for the three test cases. 
As expected the scattered field components ETSCAT and EPSCAT 
are 40 db below the source field (ETF) since in all cases R s 
is 10000. EPF is zero since the antenna radiation is 0 
polarized in this plane. Figures 5.4 through 5.6 are plots of 
field magnitude versus angle 0. 

The plots show the convergence of the calculated values as 
measured against the known analytic solution of a circular 
aperture in the far field. These results verify the program 
for the case where the radome is nearly transparent for large 
real surface resistance in accordance with Equation (1.1). 
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Figure 5.4 Electric Field Components Calculated by LDBORMM.F 
for Test 1. 
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When the radome is transparent the calculated radiation 
pattern becomes the closed form solution of the far field. 

B. EFFECTS OF RADOME ON TRANSMISSION OF ELECTRIC FIELD 

Test cases were conducted for various radome shapes and 
surface impedances. In Figures 5.7 through 5.9 the calculated 
radiation patterns with the radome are compared to the 
patterns of the isolated antenna. A comparison of the two 
curves illustrates the effect of the radome on the antenna's 
performance. The following parameters were constant for the 
three test cases: 

* resistance = 0+jl700 n; (reflection coefficient r~0.2.). 

* number of azimuthal modes = 2 . 

* number of integration points;NT=2,NPHI=40,CNRHO=20,CNPHI=40. 

* antenna radius r 0 = 2 A.. 

* radius of the BOR base was 3A (open base-no body structure) . 
The radome shapes tested were: 

* TEST 1. a sphere of radius, R=3A.. 

* TEST 2. a cone with che half angle = 30°. 

* TEST 3. an ogive with a parent circle radius R = 7.5^. 

The forward direction is the region 0 < 90° while the rear 
hemisphere is 0 > 90°. The plotted data in Figures 5.8 and 5.9 
have significant backscatter (0>9O°) due to the shape of the 
radome. The electric field magnitude of the sphere (Figure 
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5.7) for 0>9O° was below a -20 db threshold. This threshold 
was arbitr rily established for purposes of plotting the data. 
The backscatter from the sphere is negligible since the 
radiation from the antenna at each point on the surface is 
near normal incidence. The cone and the ogive shapes reflect 
energy towards other parts of the surface. These multiple 
reflections raise the sidelobes, especially at wide angles 
from the antenna main beam. 

The coupling of the reflected energy to the antenna is 
assumed to be negligible in this model. This assumption is not 
valid for geometries or dielectric materials which cause 
significant backscatter, because the reflections from the 
radome perturb the antenna excitation. Not only does this 
effect increase the sidelobes, but also causes an increase in 
the input VSWR. For a well designed radome, this mutual 
interaction is a secondary effect. 

C. RECOMMENDATIONS FOR FURTHER STUDY 

The following suggestions are made for further development 
and improvement of the program LDBORMM.F : 

1. Obtain test data for a real system. Run the program for the 
data and compare the results of the computer model to the 
actual measured test data. 

2. Develop test data from an experimental model constructed in 
the ECE Departmental Transient Electromagnetic Scattering Lab. 
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Compare the measured fields with the predictions of the 
computer model. 

3. Perform a parametric design study to access the impact of 
radome shape, source location and dielectric profiles on the 
radiation pattern. 

4. Perform an exhaustive study of the convergence properties 
of the MM solution. 

5. To improve the execution time some of the Gaussian 
integration loops can be replaced with a delta function 
approximation for the integral. 

6. Symmetry between the positive and negative modes can also 
be exploited to reduce run time. 

This thesis has developed flexible, modular code for 
modeling the electric field radiation pattern for a source 
enclosed in a dielectric body of revolution. Further work 
should concentrate on validation and improvements to reduce 
execution time. After these enhancements are incorporated, the 
code can be applied to the radome design process. 
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PHERICAL RADOME (*), AHALYT 



Figure 5.7 Effects of Radome on Far Field Radiation 
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APPENDIX A. SOURCE CODE AND DATA FILE 


A. SOURCE CODE FOR MAIN PROGRAM 


C PROGRAM : ldbormm.f 

C DATE : 23 January 1992 

C REVISED : 21 February 1992 

C PROGRAMMERS : D. JENN, R. FRANCIS 
C 

C » SPECIALIZED FOR ARBITRARY CIRCULAR APERTURE EXCITATION « 

C » EXCITATION IS SPECIFIED IN THE SUBROUTINE GENEX (....) « 
C 

c BASED ON MAUTZ AND HARRINGTON'S COMPUTER PROG FOR BORS. 
c all or part of the surface may be covered with a surface 
c impedance. 

C SOURCE IS AT BOR COORDINATE SYSTEM ORIGIN. 

C 

c imp=0 perfect conductor 
c iprint=0 print pattern points to unit 8 
c iseg=0 print the generating curve points to unit 8 


C BEGIN MAIN PROGRAM************ 

CHARACTER*8 GNT,GNPHI,CGR,CGP 
CHARACTER*14 TPTS,PPTS,PHIPTS,RHOPTS 
COMPLEX EP,ET,Z(100000),R(1600),B(1600),C(800),U,UC 
COMPLEXET1,EP1,ET2,EP2,EC,EX,ZLO(400),ZL(2400),ETF,EPF 
COMPLEX EXP1,EXP2,CONJG,CEXP,CMPLX,CT(3000),IMP,JX 
COMPLEX CIRCTHETA,CIRCPHI 

DIMENSION RH(400),ZH(400),XT(4),AT(4),IPS(800) 
DIMENSION A(100),X(100),EXP(500),ANG(500),ECP(500) 
DIMENSION XP(100),AP(100),XR(50),AR(50) 

DIMENSION ECV(500),EXV(500),PHC(500),PHX(500) 

REAL ETSCAT(500),EPSCAT(500),ETHF(500),EPHF(500) 
INTEGER NT,NPHI,CNRHO,CNPHI,NP,SELECTION 
REAL MODE,BASE,RS,ZP,RHB,ZHB 
DATA PI,START,STOP/3.1415926,0.,90./ 

DATA IPRINT/0/ 

Rad=PI/180. 

ECX=0. 

BK=2.*PI 
U=(0.,1.) 
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U0=(0.,0.) 

UC=-U/4./PI 

JK=(0.0,6.283185307) 

C Select the subroutine to generate BOR. 

WRITE(6,*)'MAKE BOR SELECTION' 

WRITE(6,*)' ' 

WRITE(6,*)'Enter number from menu to make selection' 
WRITE(6,*)' ' 

WRITE(6,*)' Selection Number BOR Geometry' 

WRITE(6,*)' 1 OGIVE' 

WRITE(6,*)' 2 SPHERE' 

READ(5,*)SELECTION 

IF(SELECTION.EQ.1)THEN 

CALL OGIVE(NP,ZH,RH,BASE,RS,ZP) 

ELSEIF(SELECTION.EQ.2)THEN 

CALL TESTSPHERE(NP,ZH,RH,BASE,RS,ZP) 

ENDIF 


C*********CALL OGIVE or TESTSPHERE. 

IF (NP.GT.399)THEN 

WRITE(6,*)'MAXIMUM NUMBER OF POINTS(NP) IS 399' 
GOTO 999 
ENDIF 

DT=90.0/FLOAT(NP-1) 

WRITE(6,*)'ENTER THE FILENAMES gaus###' 

WRITE (6, *) '.' 

WRITE(6,*)'ENTER THE FILENAME IN T (NT)' 

READ(5,*)GNT 

WRITE(6,*)'ENTER THE FILENAME IN PHI (NPHI)' 

READ(5,*)GNPHI 

WRITE(6,*)'ENTER THE FILENAME IN CNRHO' 

READ(5,*)CGR 

WRITE(6,*)'ENTER THE FILENAME IN CNPHI' 

READ(5,*)CGP 

C OPEN THE FILES FOR THE gaus/gaus### 

TPTS='gaus/'//GNT 
PPTS='gaus/'//GNPHI 
RHOPTS='gaus/'//CGR 
PHIPTS='gaus/'//CG? 

OPEN (1, FILE=TPi‘S , STATUS=' OLD ' ) 

OPEN(2,FILE=PPTS,STATUS='OLD') 

OPEN(3,FILE=RHOPTS,STATUS^'OLD') 

OPEN(4,FILE=PHIPTS,STATUS='OLD') 

READ(1,*)NT 

IF(NT.GT.4)THEN 

WRITE(6,*)'MAXIMUM NUMBER OF POINTS(NT) IS 4' 
GOTO 999 
ENDIF 

READ(2,*)NPHI 

IF(NPHI.GT.200)THEN 


70 







o o o o 


WRITE(6,*)'MAXIMUM NUMBER OF POINTS(NPHI) IS 200' 
GOTO 999 
ENDIF 

READ(3,*)CNRHO 

IF(CNRHO.GT.50)THEN 

WRITE(6,*)'MAXIMUM NUMBER OF POINTS(CNRHO) IS 50' 
GOTO 999 
ENDIF 

READ(4,*)CNPHI 

IF(CNPHI.GT.100)THEN 

WRITE(6,*)'MAXIMUM NUMBER OF POINTS(CNPHI) IS 100' 
GOTO 999 
ENDIF 

C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 

DO 1 M=1,NT 

READ(1,*,END=1)XT(M),AT(M) 

1 CONTINUE 

DO 2 M=1,NPHI 

READ(2,*,END=2) X (M) , A(M) 

2 CONTINUE 

DO 3 M=l,CNRHO 

READ(3 , * ,END=3)XR(M),AR(M) 

3 CONTINUE 

DO 4 M=l,CNPHI 

READ(4 , * ,END=4)XP(M),AP(M) 

4 CONTINUE 
CLOSE(1) 

CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 


C 


WRITE(6,*)'ENTER 
READ(5,*)MODE 
WRITE(6,*)'ENTER 
READ(5,*)P 
PHI=P*RAD 
WRITE(6,*)'ENTER 
READ(5,*)SA 
SCAN=-SA*RAD 
WRITE(6,*)'ENTER 
READ (5, *)IMP 
WRITE(6,*)'ENTER 
READ(5,*)ARAD 


MODE (FLOAT)' 

PHI (observation) IN DEGREES' 

THE SCAN ANGLE IN DEGREES' 

COMPLEX IMPEDANCE' 

THE ANTENNA RADIUS (wavelengths)' 


ENTER THE NUMBER OF AZIMUTHAL MODES 
(n=-MODES,...,0,...,+MODES) 

MODES=MODE 


71 



o o o o 


NBL0CK=2*MODES+1 
MHI=M0DES+1 

OPEN(8,FILE= 7 out1dbor') 

WRITE(8,8000) 2.*BASE,NP,RS,ZP 
8000 FORMAT(//,5X,' *** BOR RADIATION PATTERN FOR A CIRCULAR' 
*' DISC RADIATOR USING GENEX ***', 

*//,2X,'BOR DIAMETER (WAVELENGTHS)=',F5.2,/,2X, 

* 'NUMBER OF GENERATING POINTS (NP)=',I4, 

*//,2X ,'SURFACE RADIUS',F5.2,2X,'ZPRIME',F5.2) 

WRITE(8,30) NT,NPHI 

30 FORMAT(/,12X,' NT NPHI',/,10X,I5,2X,I5) 

IF(ISEG.EQ.O) WRITE(8,1300) 

1300 FORMAT(/, 10X, 'INDEX',8X, 'Z(I) ',10X, 'RHO(I) ',12X, 'SURF 

IMPED') 

DO 52 1=1,NP 

IF(ABS(ZH(I)).LT..001) ZH(I)=0. 

IF(ABS(RH(I)).LT..001) RH(I)=0. 

ZHB=ZH(I)/BK 
RHB=RH(I)/BK 
ZlO(i)=IMP 

ASSIGN SURFACE IMPEDANCE AT THIS POINT. THE SURFACE 
IMPEDANCE OF SEGMENT I IS ZLO(I) 

WRITE(8,8004) I,ZHB,RHB,ZLO(I) 

52 CONTINUE 

8004 FORMAT(11X,14,4X,F8.3,8X,F8.3,6X,2F8.2) 

C 

C MODE LOOP TO CALCULATE THE CURRENT COEFFICIENTS. POS AND NEG 
C MODES DONE IN THE SAME ITERATION OF THE LOOP 


C *************ZLOAD,ZMAT,GENEX,DECOMP,SOLVE 

IF(CABS(IMP).NE.0) CALL ZLOAD(NP,RH,ZH,ZLO,ZL) 

DO 400 M=1,MHI 

NM=M-1 

CALL ZMAT(NM,NM,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 

IF(CABS(IMP).NE.O) CALL ZTOT(MT,MP,ZL,Z) 

CALL GENEX(NM,NP,NT,NPHI,CNRHO,CNPHI,XT,AT,X,A, 
* XR,AR,XP,AP,SCAN,PHI,ARAD,RH,ZH,B) 

CALL DECOMP(N,IPS,Z) 

CALL SOLVE(N,IPS,Z,B,C) 


£★***********★**★*****★★**★■**★■*****★** 

C STORE CURRENT COEFFICIENTS IN ONE LONG COLUMN VECTOR 
NTOPl=MODES-NM 
NTOP2=NBLOCK-(NTOP1+1) 

NS2=NT0P1*N 
NS1=NT0P2*N 
C POSITIVE MODE 
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DO 401 L=1,N 
401 CT(L+NS1)=C(L) 

IF(NM.NE.O) THEN 
C NEGATIVE MODE 
NMN=-NM 


C**********************ZMAT,GENEX,DECOMP,SOLVE. 

CALL ZMAT(NMN,NMN,NP,NPHI,NT,RH, ZH,X, A,XT,AT,Z) 

IF(CABS(IMP).NE.0) CALL ZTOT(MT,MP,ZL,Z) 

CALL GENEX(NMN,NP,NT,NPHI,CNRHO,CNPHI,XT,AT , X,A, 

* XR,AR,XP,AP,SCAN,PHI,ARAD,RH,ZH,B) 

CALL DECOMP(N,IPS,Z) 

CALL SOLVE(N,IPS,Z,B,C) 
q********************************** 

DO 402 L=1,N 
402 CT(L+NS2)=C(L) 

ENDIF 

400 CONTINUE 
IT=NP 

DT=STOP/FLOAT(NP-1.) 

BEGIN PATTERN LOOP FROM START TO STOP IN INCREMENTS OF DT 
(ALL IN DEG) 

DO 500 1=1,IT 
THETA=FLOAT(1-1)*DT+START 
THX=THETA*RAD 
PHR=PHR0 
RHB=RH(I)/BK 
ZHB=ZH(I)/BK 
IF(THETA.GT.180.) THEN 
PHR=PHR0+PI 
THX=(360.-THETA)*RAD 
ENDIF 

ET1=(0.,0.) 

EP1=(0.,0.) 

ET2=(0.,0.) 

EP2=(0.,0.) 

DO 300 M=1,MHI 
NM=M-1 

EXP1=CEXP(CMPLX(0.,FLOAT(NM)*PHR)) 

EXP2=CONJG(EXP1) 


c*****************plane 

CALL PLANE(NM,NM,NP,NT,RH,ZH,XT,AT,THX,R) 

NTOPl=MODES-NM 

NTOP2=NBLOCK-(NTOP1+1) 

NS2=NT0P1*N 
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NS1=NT0P2*N 
DO 250 L=1,MT 

ET1=ET1+R(L)*CT(L+NS1)*EXP1 
EP1=EP1+R(L+N)*CT(L+NS1)*EXP1 
IF(NM.EQ.O) GO TO 250 
ET1=ET1+R(L)*CT(L+NS2)*EXP2 
EP1=EP1-R(L+N)*CT(L+NS2)*EXP2 
250 CONTINUE 

DO 260 L=1,MP 

ET2=ET2+R(L+MT)*CT(L+NS1+MT)*EXP1 
EP2=EP2+R(L+MT+N)*CT(L+NS1+MT)*EXP1 
IF(NM.EQ.O) GO TO 260 
ET2=ET2-R(L+MT)*CT(L+NS2+MT)*EXP2 
EP2=EP2+R(L+MT+N)*CT(L+NS2+MT)*EXP2 
260 CONTINUE 
300 CONTINUE 

C DISK CONTRIBUTION IN THE FAR FIELD IS ETF,EPF 
RRR=SQRT(RHB**2+ZHB**2) 

C*************CIRCTHETA,CIRCPHI 

ETF=CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 

* SCAN,PHI,RHB,ZHB) 

EPF=CIRCPHI(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 

* SCAN,PHI,RHB,ZHB) 

c i deleted the rrr*cexp(jk*rrr) factor from ETF & EPF. 

C************************************* 

ETHF(I)=CABS(ETF) 

EPHF(I)=CABS(EPF) 

C TOTAL E-THETA AND E-PHI COMPONENTS 
ET=ET3 +ET2+ETF 
EP=EP1+EP2+EPF 
ETSCAT(I)=CABS(ET-ETF) 

EPSCAT(I)=CABS(EP-EPF) 

EC=ET 

EX=EP 

ECV(I)=CABS(EC) 

EXV(I)=CABS(EX) 

ECR=REAL(EC) 

ECI=AIMAG(EC) 

EXR=REAL(EX) 

EXI=AIMAG(EX) 

PHC(I)=ATAN2(ECI,ECR+1.e-10)/RAD 
PHX(I)=ATAN2(EXI,EXR+1.e-10)/RAD 
ANG(I)=THETA 

ECX=AMAX1(ECX,ECV(I),EXV(I)) 

500 CONTINUE 

WRITE(6,*) 'MAX E VALUE= / ,ECX 
WRITE(8,103) P,ECX 

103 FORMAT(/,1OX,'PHI OF RECEIVER (DEG)=',F8.2,/,10X, 

* 'MAXIMUM FIELD VALUE (V/M)=',E15.5) 
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DO 600 1=1,IT 

ECP(I)=(ECV(I)/ECX)**2 

EXP(I)=(EXV(I)/ECX)**2 

ECP(I)=AMAX1(ECP(I),1.E-10) 

EXP(I)=AMAX1(EXP(I),1.E-10) 

ECP(I)=10.*ALOG10(ECP(I)) 

EXP(I)=10.*ALOG10(EXP(I)) 

600 CONTINUE 

IF(IPRINT.EQ.0) THEN 
WRITE(8,5015) 

5015 FORMAT(//,7X,'ANGLE',15X,'CO-POL',25X,'X-POL',/,7X, 

1'(DEG)',4X,2('(VOLTS)',4X,'(DEG)',3X,'(DB-REL)',4X)) 

DO 9000 L=1,IT 

WRITE(8,5016)ANG(L),ECV(L),PHC(L),ECP(L),EXV(L),PHX(L) 
1,EXP(L) 

5016 FORMAT(5X,F6.2,3X,2(F8.4,3X,F7.2,3X,F7.2,3X)) 

9000 CONTINUE 

ENDIF 

OPEN(2,file='bang.m') 

OPEN (3 , f ile='bcpole .in') 

OPEN(4,file='bxpole.m') 

OPEN(7,file='etscat.m') 

OPEN(8,file='epscat.m') 

OPEN(9,file='etf.m') 

OPEN(10,file='epf.m') 

DO 9097 1=1,IT 
WRITE(2,5019) ANG(I) 

WRITE(3,5019) ECV(I) 

WRITE(4,5019) EXV(I) 

C change ecv,ecv to ecp & exp to get normalized db values. 
WRITE(7,5019) ETSCAT(I) 

WRITE(8,5019) EPSCAT(I) 

WRITE(9,5019) ETHF(I) 

9097 WRITE(10,5019) EPHF(I) 

5019 FORMAT(F8.3) 

CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

CLOSE(7) 

CLOSE(8) 

CLOSE(9) 

CLOSE(10) 

999 CONTINUE 
STOP 
END 

C**************END OF MAIN PROGRAM. 
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c 

C**************SUBROUTINE ZMAT. 

C REFERENCE: AN IMPROVED E-FIELD SOLUTION FOR CONDUCTING BOR 
C J.R.MAUTZ AND R.F. HARRINGTON 

C TECHNICAL REPORT TR-80-1 

C ROME AIR DEVELOPMENT CENTER 

C CONTRACT NO. F-30602-79-C-0011 

SUBROUTINE ZMAT(Ml,M2,NP,NPHI,NT,RH,ZH,X,A,XT,AT,Z) 

COMPUTE THE MM IMPEDANCE MATRIX ELEMENTS. THIS IS FROM 
MAUTZ AND HARRINGTON (NO CHANGES). 


COMPLEXZ(100000),U1,U2,U3,U4,U5,U6,U7,U8,U9,UA,UB,G4A(4),G5A 
(4) 

COMPLEXG6A(4) ,G4B(4) f G5B(4) ,G6B(4) ,H4A,H5A,H6A,H4 B,H5B 
COMPLEX H6B,UC,UD,GA(100),GB(100) 

DIMENSION RH(400),ZH(400),X(100),A(100),XT(4),AT(4) 
DIMENSION RS(400),ZS(400),D(400),DR(400),DZ(400) 
DIMENSION DM(400), C2(100),C3(100),R2(4),Z2(4) 

DIMENSION C4(100),C5(100),C6(100),Z7(4),R7(4),Z8(4),R8(4) 
CT=2 . 

CP=. 1 

DO 10 1=2,NP 
12=1-1 

RS(12) =.5 *(RH(I)+RH(12)) 

ZS(12)=.5*(ZH(I)+ZH(I2)) 

Dl=.5*(RH(I)-RH(12)) 

D2=.5*(ZH(I)-ZH(12)) 

D(12)=SQRT(D1*D1+D2*D2) 

DR(12)=D1 
DZ(12)=D2 

IF(RS(12).EQ.0.) RS(I2)=1. 

DM(12)=D(12)/RS(12) 

10 CONTINUE 
M3=M2-M1+1 
M4=M1-1 
PI2=1.570796 
DO 11 K=1,NPHI 
PH=PI2*(X(K)+: 

C2(K)=PH*PH 
SN=SIN(.5*PH 
C3(K)=4.*SN*SN 
A1=PI2*A(K) 

D4=.5*A1*C3(K) 

D5=A1*C0S(PH) 

D6=A1*SIN(PH) 

M5=K 

DO 29 M=1,M3 
PHM=(M4+M)* PH 
A2=COS(PHM) 
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C4(M5)=D4*A2 
C5(M5)=D5*A2 
C6(M5)=D6*SIN(PHM) 
M5=M5+NPHI 
29 CONTINUE 

11 CONTINUE 
MP=NP-1 
MT=MP-1 
N=MT+MP 
N2N=MT*N 
N2=N*N 
U1=(0.,.5) 

U2=(0.,2.) 

JN=-1-N 

DO 15 JQ=1,MP 
KQ=2 

IF(JQ.EQ.l) KQ=1 
IF(JQ.EQ.MP) KQ=3 
R1=RS(JQ) 

Z1=ZS(JQ) 

D1=D(JQ) 

D2=DR(JQ) 

D3=DZ(JQ) 

D4=D2/R1 
D5=DM(JQ) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6+D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO 12 L=1,NT 

R2(L)=R1+D2*XT(L) 

Z2(L)=Z1+D3*XT(L) 

12 CONTINUE 
U3=D2*U1 
U4=D3*U1 

DO 16 IP=1,MP 
R3=RS(IP) 

Z3=ZS(IP) 

R4=R1-R3 
Z4=Z1-Z 3 
FM=R4*SV+Z4*CV 
PHM=ABS(FM) 

PH=ABS(R4*CV-Z4*SV) 
D6=PH 

IF(PHM.LE.Dl) GO TO 26 
D6=PHM-D1 

D6=SQRT(D6*D6+PH*PH) 

26 IF(IP.EQ.JQ) GO TO 27 
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KP=1 

IF(T6.GT.D6) KP=2 
IF(R6.GT.D6) KP=3 
GO TO 28 

27 KP=4 

28 GO TO (41,42,41,42),KP 
42 DO 40 L=1,NT 

D7=R2(L)-R3 
D8=Z2(L)-Z3 
Z7(L)=D7*D7+D8*D8 
R7(L)=R3*R2(L) 

Z8(L) = .25*Z7(L) 

R8(L)=.25*R7(L) 

40 CONTINUE 

Z4=R4*R4+Z4 *Z4 
R4=R3*R1 
R5=.5*R3*SV 
DO 33 K=1,NPHI 
A1=C3(K) 

RR=Z4+R4*A1 
UA=0. 

UB=0. 

IF(RR.LT.T62) GO TO 34 
DO 35 L=1,NT 
R=SQRT(Z7(L)+R7(L)*A1) 

SN=-SIN(R) 

CS=COS(R) 

UC=AT(L)/R*CMPLX(CS,SN) 
UA=UA+UC 
UB=XT(L)*UC+UB 
35 CONTINUE 
GO TO 36 
34 DO 37 L=1,NT 

R=SQRT(Z8(L)+R8(L)*A1) 
SN=-SIN(R) 

CS=COS(R) 

UC=AT(L)/R*SN*CMPLX(-SN,CS) 

UA=UA+UC 

UB=XT(L)*UC+UB 

37 CONTINUE 
A2=FM+R5*A1 
D9=RR-A2 *A2 
R=ABS(A2) 

D7=R-D1 

D8=R+D1 

D6=SQRT(D8*D8+D9) 

R=SQRT(D7*D7+D9) 

IF(D7.GE.0.) GO TO 38 
Al=ALOG((D8+D6)*(-D7+R)/D9)/D1 
GO TO 39 

38 Al=ALOG((D8+D6)/(D7+R))/Dl 
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39 UA=A1+UA 

UB=A2 *(4./(D6+R)-A1)/Dl+UB 
36 GA(K)=UA 
GB(K)=UB 
33 CONTINUE 
K1=0 

DO 45 M=1,M3 
H4A=0. 

H5A=0. 

H6A=0. 

H4B=0. 

H5B=0. 

H6B=0. 

DO 46 K=1,NPHI 
K1=K1+1 
D6=C4(Kl) 

D7=C5(K1) 

D8=C6(Kl) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
H5B=D7 *UB+H5B 
H6B=D8*UB+H6B 
46 CONTINUE 
G4A(M)=H4A 
G5A(M)«Ii5A 
G6A(M)=H6A 
G4B(M)=H4B 
G5B(M)=H5B 
G6B(M)=H6B 
45 CONTINUE 

IF(KP.NE.4) GO TO 47 
A2=D1/(PI2*R1) 

D6=2./D1 
D8=0. 

DO 63 K=1,NPHI 
A1=R4 *C2(K) 

R=R4 *C3(K) 

IF(R.LT.T62) GO TO 64 
D7=0. 

DO 65 L=1,NT 

D7=D7+AT(L)/SQRT(Z7(L)+A1) 

65 CONTINUE 
GO TO 66 

64 A1=A2/(X(K)+1.) 

D7=D6*ALOG(Al+SQRT(1.+A1*A1)) 

66 D8=D8+A(K)*D7 
63 CONTINUE 
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A1=.5*A2 

A2=1./A1 

D8=-PI2*D8+2./Rl*(BLOG(A2)+A2*BLOG(Al)) 

DO 67 M=1,M3 
G5A(M)=D8+G5A(M) 

67 CONTINUE 
GO TO 47 
41 DO 25 M=1,M3 
G4A(M)=0. 

G5A(M)=0. 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)=0. 

G6B(M)=0. 

25 CONTINUE 

DO 13 L=1,NT 
A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO 17 K=1,NPHI 

R=SQRT(Z4+R4*C3(K)) 

SN=-SIN(R) 

CS=COS(R) 

GA(K)=CMPLX(CS f SN)/R 
17 CONTINUE 
D6=0. 

IF(R62.LE.Z4) GO TO 51 

DO 62 K=1,NPHI 

D6=D6+A(K)/SQRT(Z4+R4 *C2(K)) 

62 CONTINUE 

Z4=3.141593/SQRT(Z4/R4) 

D6=-PI2*D6+ALOG(Z4+SQRT(l.+Z4*Z4))/SQRT(R4) 
51 A1=AT(L) 

A2=XT(L)*A1 
K1=0 

DO 30 M=1,M3 
U5=0. 

U6=0. 

U7=0. 

DO 32 K=1,NPHI 
UA=GA(K) 

K1=K1 4 1 

U5=C4(Kl)*UA+U5 
U6=C5(Kl)*UA+U6 
U7=C6(Kl)*UA+U7 
32 CONTINUE 
U6=D6+U6 

G4A(M)=A1*U5+G4A(M) 

G5A(M)=A1*U6+G5A(M) 
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G6A(M)=A1*U7+G6A(M) 
G4B(M)=A2*U5+G4B(M) 
G5B(M)=A2*U6+G5B(M) 
G6B(M)=A2*U7+G6B(M) 

3 0 CONTINUE 

13 CONTINUE 

47 A1=DR(IP) 

UA=A1*U3 
UB=DZ(IP)*U4 
A2=D(IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JM=JN 

DO 31 M=1,M3 
FM=M4+M 
A1=FM*DM(IP) 

H5A=G5A(M) 

H5B=G5B(M) 

H4A=G4A(M)+H5A 
H4B=G4B(M)+H5B 
H6A=G6A(M) 

H6B=G6B(M) 

U7=UA*H5A+UB*H4A 

U8=UA*H5B+UB*H4B 

U5=U7-U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7*(H6A+D4 *H6B) 

UD=FM*D5*H4A 

K1=IP+JM 

K2=K1+1 

K3=K1+N 

K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO TO (18,20,19),KQ 

18 Z(K6)=U8+U9 
IF(IP.EQ.l) GO TO 21 
Z(K3)=Z(K3)+U6-U7 

Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 

21 Z(K4)=U6+U7 
Z(K8)=UC+UD 
GO TO 22 

19 Z(K5)=Z(K5)+U8-U9 
IF(IP.EQ.l) GO TO 23 
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non 


Z (K1)=Z (Kl) +U5+L 7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 

23 Z (K2 )=Z(K2)+U5-U7 
Z(K8)=UC+UD 

GO TO 22 

20 Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.l) GO TO 24 
Z(K1)=Z(Kl)+U5+U7 
Z(K3)=Z(K3)+U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 

24 Z(K2)=Z(K2)+U5-U7 
Z(K4)=U6+U7 
Z(K8)=UC+UD 

22 Z(K8+MT)=U2*(D8*(H5A+D4*H5B)-A1*UD) 
JM=JM+N2 

31 CONTINUE 

16 CONTINUE 
JN=JN+N 

15 CONTINUE 
RETURN 
END 


C************* **SUBROUTINE SOLVE. 

C REFERENCE : TECHNICAL REPORT TR-80-1 
SUBROUTINE SOLVE(N,IPS,UL,B,X) 

SEE MAUTZ & HARRINGTON FOR DETAILS 

COMPLEX UL( 100000) , B ( 800) , X (800) . c ’ r ' v * 
DIMENSION IPS(800) 

NP=N+1 
IP=IPS(1) 

X (1) =B (IP) 

DO 2 1=2,N 
IP=IPS(I) 

IPB=IP 
IM1=I-1 
SUM=(0.,0.) 

DO 1 J=1,IM1 
SUM=SUM+UL(IP)*X(J) 

1 IP=IP+N 

2 X(I)=B(IPB)-SUM 
K2=N*(N-l) 

IP=IPS(N)+K2 
X(N)=X(N)/UL(IP) 

DO 4 IBACK=2,N 
I=NP-IBACK 


82 




non 


K2=K2-N 
IPI=IPS(I)+K2 
IP1=I+1 
SUM=(0.,0.) 

IP=IPI 

DO 3 J=IP1,N 
IP=IP+N 

3 SUM=SUM+UL(IP)*X(J) 

4 X(I)=(X(I)-SUM)/UL(IPI) 
RETURN 

END 


C* * * * * *SUBROUTINE DECOMP 
C REFERENCE: TECHNICAL REPORT TR-80-1 
SUBROUTINE DECOMP(N,IPS,UL) 

SEE MAUTZ & HARRINGTON FOR DETAILS 


COMPLEX UL(IOOOOO),PIVOT,EM 
DIMENSION SCL(800),IPS(800) 

DO 5 1=1,N 
IPS(I)=I 
RN=0. 

J1=I 

DO 2 J=1,N 

ULM=ABS(R£AL(UL(J1)))+ABS(AIMAG(UL(Jl))) 

J1=J1+N 

IF(RN-ULM) 1,2,2 

1 RN=ULM 

2 CONTINUE 
SCL(I)=1./RN 

5 CONTINUE 
NM1=N-1 
K2=0 

DO 17 K=1,NM1 
BIG=0. 

DO 11 I=K,N 
IP=IPS(I) 

IPK=IP+K2 

SIZE=(ABS(REAL(UL(IPK)))+ABS(AIMAG(UL(IPK))))*SCL(IP) 
IF(SIZE-BIG) 11,11,10 

10 BIG=SIZE 
IPV=I 

11 CONTINUE 

IF(IPV-K) 14,15,14 

14 J=IPS(K) 

IPS(K)=IPS(IPV) 

IPS(IPV)=J 

15 KPP=IPS(K)+K2 
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PIVOT=UL(KPP) 

KP1=K+1 
DO 16 I=KP1,N 
KP=KPP 
IP=IPS(I)+K2 
EM=-UL(IP)/PIVOT 
18 UL(IP)=-EM 
DO 16 J=KP1,N 
IP—IP+N 
KP=KP+N 

UL(IP)=UL(IP)+EM*UL(KP) 

16 CONTINUE 
K2=K2+N 

17 CONTINUE 
RETURN 
END 


C************FUNCTION BLOG 

C REFERENCE: TECHNICAL REPORT TR-80-1 ;(page 56) 
FUNCTION BLOG(X) 

IF(X.GT..1) GO TO 1 
X2=X*X 

BLOG=((.075*X2-.1666667)*X2+1.)*X 
RETURN 

1 BLOG=ALOG(X+SQRT(1.+X*X)) 

RETURN 

END 


C**************** * *SUBROUTINE PLANE 

C REFERENCE: TECHNICAL REPORT TR-80-1 ;(pages 57-64) 
SUBROUTINE PLANE(Ml,M2,NP,NT,RH,ZH,XT,AT,THR,R) 

PLANE WAVE EXCITATION VECTOR. FROM MAUTZ AND HARRINGTON. 
MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 

COMPLEX R(1600),U,U1,UA,UB,FA(5),FB(5),F2A,F2B,FIA,FIB 
COMPLEX U2,U3,U4,U5 

DIMENSION RH(400),ZH(400),XT(4),AT(4),R2(4) 

DIMENSION Z2(4),BJ(2000) 

MP=NP-1 
MT=MP-1 
N=MT+MP 
N2=2*N 
CC=COS(THR) 

SS=SIN(THR) 

U=(0.,1.) 

U1=3.141593*U**M1 
M3=M1+1 
M4=M2 +3 
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IF(Ml.EQ.O) M3=2 

M5=Ml+2 

M6=M2+2 

DO 12 IP=1,MP 

K2=IP 

I»IP+1 

DR=.5 *(RH(I)-RH(IP)) 

DZ=.5*(ZH(I)-ZH(IP)) 

D1=SQRT(DR*DR+DZ*DZ) 

Rl=.25*(RH(I)+RH(IP)) 

IF(R1.EQ.0.) Rl=l. 

Zl=.5*(ZH(I)+ZH(IP)) 

DR=.5*DR 

D2=DR/R1 

DO 13 L=1,NT 

R2(L)=R1+DR*XT(L) 

Z2(L)=Z1+DZ*XT(L) 

13 CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO 23 M=M3 , M4 
FA(M)=0. 

FB(M)=0. 

23 CONTINUE 

DO 15 L=1,NT 
X=SS*R2(L) 

IF(X.GT..5E-7) GO TO 19 
DO 20 M=M3,M4 
BJ(M)=0. 

20 CONTINUE 
BJ (2) =1. 

S=1. 

GO TO 18 

19 M=2.8*X+14--2./X 

IF(X.LT..5) M=ll.8+ALOG10(X) 
IF(M.LT.M4) M=M4 
BJ (M) =0. 

JM=M-1 
BJ(JM)=1. 

DO 16 J=4,M 
J2=JM 
JM=JM-1 
J1=JM-1 

BJ(JM)=J1/X*BJ(J2)-BJ(JM+2) 

16 CONTINUE 
S=0. 

IF(M.LE.4) GO TO 24 
DO 17 J=4,M,2 
S=S+BJ(J) 

17 CONTINUE 
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24 S=BJ(2)+2.*S 
18 ARG=Z2(L)*CC 

UA=AT(L)/S*CMPLX(COS(ARG),SIN(ARG)) 
UB=XT(L)*UA 
DO 25 M=M3,M4 
FA (M) =BJ (M) *UA+FA (M) 

FB(M)=BJ(M)*UB+FB(M) 

25 CONTINUE 
15 CONTINUE 

IF(Ml.NE.O) GO TO 26 
FA(1)=FA(3) 

FB(1)=FB(3) 

26 UA=U1 

DO 27 M=M5,M6 

M7=M-1 

M8=M+1 

F2A=UA*(FA(M8)+FA(M7)) 

F2B=UA*(FB (M8)+FB(M7)) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7)) 
F1B=UB*(FB(M8)-FB(M7)) 

U4=D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+MT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(K5)=U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 
12 CONTINUE 

RETURN 

END 


C*************SUBR0UTINE ZLOAD. 

C REFERENCE: COMPUTATION OF RADIATION AND SCATTERING FOR 
C LOADED BODIES OF REVOLUTION 

C HARRINGTON AND MAUTZ 

C REPORT AFCRL-70-0046 

C AIRFORCE CAMBRIDGE RESEARCH LABS 

C CONTRACT NO. F-19628-68-C-0180 
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SUBROUTINE ZLOAD(NP,RH,ZH,ZO,Z) 


C 

C COMPUTES IMPEDANCE MATRIX ELEMENTS FOR LOADED BODIES OF REV 
C Z0(I) IS THE SURF IMPEDANCE OF THE ITH SEGMENT (NP-1 
C SEGMENTS) Z(.) ARE THE IMPEDANCE MATRIX TERMS (TRIDIAGONAL 
C FOR T-T SUBMATRIX; DIAGONAL FOR P-P SUBMATRIX). STORED IN 
C COL VECTOR. 

C 


COMPLEXC1,C2,ZO(400), Z(2400),XI,X2,X3,Y1,Y2,Y3,FN(400) 
COMPLEX U1,U2,U3,XI,YI 

DIMENSION RH(400) f ZH(400),RS(400),D(400),SV(400) 

PI=3.14159 

MT=NP-2 

MP=NP-1 

N=MT+MP 

DO 10 IP=2,NP 

II=IP-1 

DR=RH(IP)-RH(II) 

DZ=ZH(IP)-ZH(II) 

D(II)=SQRT(DR*DR+DZ *DZ) 

SV(II)=DR/D(II) 

RS(II)=.5*(RH(IP)+RH(II)) 

DS=D(II)*SV(II)/2. 

Q1=RS(II)+DS 
Q2=RS(II)-DS 
FN(II)=1. 

IF((ABS(Q2).GT.l.E-6).AND.(ABS(Ql).GT.l.E-6)) 

* FN(II)=ALOG(Q1/Q2) 

10 CONTINUE 
L0=MT*3-2 
DO 20 1=1,MP 
C1=PI*Z0(I) 

IF(I.EQ.MP) GO TO 80 
KI=2 

IF(I.EQ.l) KI=1 
IF(I.EQ.MT) KI=3 
11=1+1 
C2=PI*Z0(II) 

A=SV(I) 

IF(ABS(A).LT.l.E-6) GO TO 41 
Xl=Cl*FN(I)/2./A 

X2=C1*2./A*(1.-RS(I)*FN(I)/D(I)/A) 

X3=-X2*RS(I)/D(I)/A 
GO TO 42 

41 CONTINUE 
Xl=Cl/2./RS(I)*D(I) 

X2=(0.,0.) 

X3=Cl*D(I)/6./RS(I) 

42 CONTINUE 
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A=SV(II) 

IF(ABS(A).LT.l.E-6) GO TO 45 
Yl=C2*FN(II)/2./A 

Y2=C2*2./A*(1.-RS(II)*FN(II)/D(II)/A) 
Y3=-Y2*RS(II)/D(II)/A 
GO TO 40 
45 CONTINUE 

Yl=C2/2./RS(II)*D(II) 

Y2=(0.,0.) 

Y3=C2*D(II)/6./RS(II) 

40 CONTINUE 
C 

C DEFINE TRIDIAGONAL ELEMENTS FOR T-T SUBMATRIX (STORED IN 
C COLS) (Ul- DIAG; U2- LOWER; U3- UPPER) 

C 

XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(KI.EQ.l) XI=C1/SV(I) 

C IF(KI.EQ.3) YI=C2/SV(II) 

U1=XI+YI 

U2=X1-X3 

U3=Y1-Y3 

L=2+(1-2)*3 

IF(KI.EQ.l) Lr=0 

L1=L+1 

L2=L+2 

L3=L+3 

go to (50,60,70),ki 
50 Z(LI)=U1 
Z(L2)=U2 
GO TO 8 0 
60 Z(L1)=U3 
Z(L2)=U1 
Z(L3)=U2 
GO TO 80 
70 Z(L1)=U3 
Z(L2)=U1 

80 Z(L0+I)=2.*C1*D(I)/RS(I) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE ZTOT(MT,MP,ZL,Z) 

C 

C ADDS THE SURF IMPEDANCE TERMS TO THE TRIDIAGONAL ELEMENTS OF 
C THE BOR IMPEDANCE MATRIX Z. 

C 

COMPLEX ZL(2400),Z(100000) 

N=MT+MP 
M0=MT*3-2 
DO 100 1=1,MP 
L0=MT*N+(1-1)*N+MT 
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IF(I.EQ.MP) GO TO 80 
KI=2 

IF(I.EQ.l) KI=1 
IF(I.EQ.MT) KI=3 
L2=(1-1)*N+I 
L1=L2-1 
L3=L2+1 
M=2+3*(1-2) 

IF(KI.EQ.l) M=0 

M1=M+1 

M2=M+2 

M3=M+3 

go to (50,60,70),ki 
50 Z(L2)=Z(L2)+ZL(M1) 
Z(L3)=Z(L3)+ZL(M2) 

GO TO 80 

60 Z(L1)=Z(L1)+ZL(M1) 
Z(L2)=Z(L2)+ZL(M2) 

Z(L3)=Z(L3)+ZL(M3) 

GO TO 80 

70 Z(LI)=Z(LI)+ZL (Ml) 
Z(L2)=Z(L2)+ZL(M2) 

80 Z(LO+I)=Z(LO+I)+ZL(MO+I) 
100 CONTINUE 
RETURN 
END 


C*********SUBROUTINE OGIVE 
C SUBROUTINE : OGIVE 
C DATE : 4 SEPTEMBER 1991 

C PROGRAMMER : R.M. FRANCIS 
C REVISED : 26 JANUARY 1992 

C COMMENTS : 

C THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION 
C (BOR) IN THE FORM OF AN OGIVE. 

C DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 

C ZH = Z CO-ORDINATE * 2*PI 
C RH = RADIUS *2*PI 

SUBROUTINE OGIVE(NP,ZH,RH,BASE) 

C NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 
C ZP = ZPRIME, THE POSITION ON Z WHERE THE RADIUS OF 
C CURVATURE STARTS 
C BASE = BASE RADIUS 

C RS = RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 
INTEGER I,NP 

REAL RH(400),ZH(400),ZP,BASE,RS,ZCOORD,RADIUS,AL 
PI=3.1415926 
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C INPUT THE VARIABLES FOR THE OGIVE,ZP,B,RS,NP 

WRITE(6,*) 'ENTER SURFACE CURVATURE (wavelengths)' 

READ(5,*) RS 

WRITE (6, *) 'ENTER ZPRIME, WHERE CURVATURE STARTS (wavelengths) ' 
READ(5,*) ZP 

WRITE(6,*) 'ENTER BASE RADIUS (wavelengths)' 

READ(5,*) BASE 

C PERFORM CALCULATIONS 

AL=SQRT(2.*BASE*RS-BASE**2) 

ZMAX= AL + ZP 
ANG=ASIN(AL/RS) 

L=ANINT((RS*ANG+ZP)*5) 

NP=2*L+1 

WRITE(6,*)'NUMBER OF POINTS IS: ',NP 
DZ= ZMAX/FLOAT(NP-1) 

DO 10 1=1,NP 

ZCOORD= (1-1)*DZ 
ZH(I)=2*PI*ZCOORD 

IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 

ELSE 

RADIUS=SQRT(RS**2-(ZCOORD-ZP)**2)+(BASE-RS) 
ENDIF 

RH(I)=2*PI*RADIUS 
10 CONTINUE 
RETURN 
END 


C*********SUBR0UTINE GENEX. 

C SUBROUT INI' : GENEX 
C DATE : 14 JANUARY 1992 

C REVISED : 17 February 1992 

C PROGRAMMER : R.M. FRANCIS 

SUBROUTINE GENEX(NM,NP,NT,NPHI,CNRHO,CNPHI,XT,AT,X,A, 

* XR,AR,XP,AP,SCAN,PHI,ARAD,RH,ZH,R) 

C SOURCE ON Z AXIS AT Z=0. 

C THE EXCITATION VECTOR IS COMPUTED FOR THE GIVEN R,THETA AND 
C PHI COMPONENTS SPECIFIED IN FUNCTION SUBROUTINES 
C CIRCTHETA,CIRCRHO,CIRCPHI. 

COMPLEX R(1600),CEXP,PSI,ST1,ST2,ST3,ST4,SP,SI,S2,S3,S4,S5 
COMPLEX CIRCTHETA,CIRCRHO,CIRCPHI 

DIMENSION RH(400),ZH(400),XT(4),AT(4),X(100),A(100) 
REAL XP(100),AP(100),XR(50),AR(50),SCAN,PHI,ARAD,RHB,ZHB 
INTEGER NM,NP,NT,NPHI,CNRHO,CNPHI,I,N,MP,MT 
PI=3.141592654 





BK=2.*PI 
MP=NP-1 
MT=MP-1 
N=MT+MP 

C LIMITS ON PHI INTEGRATION 
Pl=(2.*PI-0.)/2. 

P2=(2.*PI+0.)/2. 

DO 30 IP=1,MP 

C QUANTITIES FOR THE FIRST SEGMENT (POSITIVE SLOPE) 

I=IP+1 

II=IP 

DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D1=SQRT(DR*DR+DZ*DZ) 

Rl=.5*(RH(I)+RH(II)) 

Zl=.5*(ZH(I)+ZH(II)) 

SVP1=DR/D1 

CVP1=DZ/D1 

V1=ATAN2(DR,DZ+1.E-5) 

C QUANTITIES FOR THE SECOND SEGMENT (NEGATIVE SLOPE) 

C (SKIP IF LAST SEGMENT) 

I=IP+2 

II=IP+1 

DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D2=SQRT(DR*DR+DZ*DZ) 

R2=.5*(RH(I)+RH(II)) 

Z2=.5*(ZH(I)+ZH(II)) 

SVP2=DR/D2 

CVP2=DZ/D2 

V2=ATAN2(DR,DZ+l.E-5) 

C BEGIN PHI INTEGRATION: R TERMS IN SI AND S2; THETA TERM IN 
C S3 AND S4; PHI TERM IN S5. 

ST1=(0.,0.) 

ST2=(0.,0.) 

ST3=(0.,0.) 

ST4=(0.,0.) 

SP=(0.,0.) 

DO 20 J=1,NPHI 
PH=P1*X(J)+P2 

PSI=CEXP(CMPLX(0.,-MODE*PH)) 

S1=(0.,0.) 

S2=(0.,0.) 

S3=(0.,0.) 

S4=(0.,0.) 

S5=(0.,0.) 

IF(IP.LE.MT) THEN 

C t-CURRENT INTEGRATION FOR THE POSITIVE SLOPE 
DO 13 L=1,NT 
TP=D1*XT(L)/2. 

RHB=(R1+TP*SVP1)/BK 
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ZHB=(Z1+TP*CVP1)/BK 
TH=ATAN2(RHB,ZHB+1.E-5) 

CC=COS(Vl-TH) 

SS=SIN(Vl-TH) 

RR=SQRT(RHB**2+ZHB**2) 

S1=S1+AT(L)*(.5+TP/D1)*CC*CIRCRHO(CNPHI,XP,AP,CNRHO f XR,AR, 

* ARAD,SCAN,PHI,RHB,ZHB) 

S2=S2+AT(L)*(.5+TP/D1)*SS*CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR, 

* ARAD,SCAN,PHI,RHB,ZHB) 

13 CONTINUE 
Sl=Sl*Dl/2. 

S2=S2*Dl/2. 

C t-CURRENT INTEGRATION FOR THE NEGATIVE SLOPE 
DO 14 L= 1, NT 
TP=D2 *XT(L)/2. 

RHB=(R2+TP*SVP2)/BK 
ZHB=(Z2+TP*CVP2)/BK 
TH=ATAN2(RHB,ZHB+1.E-5) 

SS=SIN(V2-TH) 

CC=COS(V2-TH) 

RR=SQRT(RHB**2+ZHB**2) 

S 3=S 3+AT(L)*(.5-TP/D2)* CC*CIRCRHO(CNPHI,XP,AP,CNRHO,XR,AR, 

* ARAD,SCAN,PHI,RHB,ZHB) 

S 4 =S 4 +AT(L)*(.5-TP/D2)* SS * CIRCTHETA(CNPHI,XP,AP,CNRHO,XR,AR, 

* ARAD,SCAN,PHI,RHB,ZHB) 

14 CONTINUE 
S3=S3*D2/2. 

S4=S4 *D2/2. 

ENDIF 

C phi-CURRENT INTEGRATION 
DO 15 L=1,NT 
TP=D1*XT(L)/2. 

RHB=(R1+TP*SVP1)/BK 
ZHB=(Z1+TP*CVP1)/BK 
TH=ATAN2(RHB,ZHB+1.E-5) 

RR=SQRT(RHB**2+ZHB**2) 

S5=S5+AT(L)+CIRCPHI(CNPHI,XP,AP,CNRHO,XR,AR, 

* ARAD,SCAN,PHI,RHB,ZHB)/Rl*(RHB*BK) 

15 CONTINUE 
S5=S5*Dl/2. 

SP=SP+A(J)*PSI*S5 
IF(IP.LE.MT) THEN 
ST1=ST1+A(J)*PSI*S1 
ST2=ST2+A(J)*PSI*S2 
ST3=ST3+A(J)*PSI*S3 
ST4=ST4+A(J)*PSI*S4 
ENDIF 
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20 CONTINUE 

C COMPONENTS ARE STORED IN A COLUMN VECTOR: VT(1,MT), 
C VP(MT+1,N) 

IF(IP.LE.MT) THEN 

R(IP)=(ST1+ST2+ST3+ST4)*P1 

ENDIF 

R(IP+MT)=SP*P1 
30 CONTINUE 

RETURN 
END 


C * * * * * * *FUNCTION CIRCTHETA. 

C FUNCTION : CIRCTHETA 

C DATE : 1 OCTOBER 1991 

C REVISED : 17 February 1992 

C PROGRAMMER : R.M. FRANCIS 

C COMMENTS : THIS FUNCTION WILL RIGOROUSLY CALCULATE THE 
C ELECTRIC FIELD FOR A CIRCULAR APERTURE. THE FIELD IS 
C CALCULATED AT THE BOUNDARY DEFINED BY AN OGIVE. THE APERTURE 
C IS LOCATED AT Z = 0, AND IS BORE SIGHTED ON THE Z AXIS. THE 
C SUBROUTINE OGIVE IS THE SOURCE OF THE GEOMETRIC DATA 
C REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 

C ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 


COMPLEX FUNCTION CIRCTHETA (CNPHI, XP, AP, CNRHO, XR, AR, ARAD, SCAN, 

* PHI,RHB,ZHB) 

INTEGER CNPHI,CNRHO,M,N 

REAL XP(100),AP(100),XR(50),AR(50), 

* ARAD,SCAN,PHI,RHB,ZHB 
C LOCAL VARIABLES 

REAL Z,RZ,PHIPRIME,PI,S,C, 

* R,RP,X2,Yl,Y2,K 
COMPLEX J,JK,G1,G2,XI,CON,CC, 

* sumx,sumy,sumz,CIRCTHETA 
RR1=(ARAD+0.)/2. 

RR2=(ARAD-0.)/2. 

PI=3.141592654 
K=6.283185307 
J = (0.0,1.0) 

JK= (0.0,6.283185307) 

CON=(0.0,-15.0) 

CC= C0N*RR1 

C OUTER LOOP: INTEGRATE OVER PHI. . .-PKPHKPI. 

C INNER LOOP: INTEGRATE OVER RHO... 0 <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
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sumx=(0.0,0.0) 
sumy= ( 0 . 0 , 0 . 0 ) 
sumz=(0.0,0.0) 

DO 50 M=1,CNPHI 

PHIPRIME=PI*XP(M) 

DO 60 N=1,CNRHO 
RP=RR1*XR(N)+RR2 

R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
Gl=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3 *JK*R)-(K*R)**2)/R**5 
X1=JK*((RP*COS(PHIPRIME)*SIN(SCAN))-R) 

Yl=(RZ*COS(PHI))-(RP*COS(PHIPRIME)) 

X2=Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 

sumx=sumx+(CC*AP(M)*AR(N)*(G1+(X2*G2))*CEXP(XI)*RP) 
sumy=sumy+(CC*AP(M)*AR(N)*Y1*Y2*G2*CEXP(X1)*RP) 
sumz=sumz+(CC*AP(M)*AR(N)*Y1*Z*G2*CEXP(X1)*RP) 

60 CONTINUE 
50 CONTINUE 

CIRCTHETA=(C*COS(PHI)* sumx) + (C*SIN(PHI)* sumy)-(S *sumz) 

RETURN 

END 


C***********FUNCTION CIRCRHO. 

C FUNCTION : CIRCRHO 

C DATE : 1 OCTOBER 1991 

C REVISED : 17 February 1992 

C PROGRAMMER : R.M. FRANCIS 


COMPLEX FUNCTION CIRCRHO(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 

* SCAN,PHI,RHB,ZHB) 

INTEGER CNPHI,CNRHO,M,N 

REAL XP(100),AP(100),XR(50),AR(50), 

* ARAD,SCAN,PHI,RHB,ZHB 
C LOCAL VARIABLES 

REAL Z,RZ,PHIPRIME,PI,S,C, 

* R,RP,X2,Y1,Y2,K 
COMPLEX J,JK,G1,G2,X1,C0N,CC, 

* sumx,sumy,sumz,CIRCRHO 
RR1=(ARAD-0.)/2. 

RR2=(ARAD+0.)/2. 

PI=3.141592654 
K=6.283185307 
J = (0.0,1.0) 

JK= (0.0,6.283185307) 

CON=(0.0,-15.0) 

CC=C0N*RR1 

C OUTER LOOP: INTEGRATE OVER PHI. . . -PKPHKPI. 





C INNER LOOP: INTEGRATE OVER RHO... 0 <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
sumx=(0.0,0.0) 
sumy=(0.0,0.0) 
sumz=(0.0,0.0) 

DO 50 M=1,CNPHI 

PHIPRIME=PI*XP(M) 

DO 60 N=1,CNRHO 
RP=RR1*XR(N)+RR2 

R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
Gl=(((K*R)**2)-l-(JK*R))/R**3 
G2=(3+(3 *JK*R)-(K*R)* *2)/R**5 
X1=JK*((RP*COS(PHIPRIME)*SIN(SCAN))-R) 

Yl=(RZ*COS(PHI))-(RP*COS(PHIPRIME)) 

X2=Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 

sumx=sumx+(CC*AP(M)*AR(N)*(G1+(X2*G2))*CEXP(X1)*RP) 
sumy=sumy+ (CC*AP(M) *AR(N) *Y1*Y2*G2*CEXP(X1) *RP) 
sumz=sumz+(CC*AP(M)*AR(N)*Y1*Z*G2*CEXP(X1)*RP) 

60 CONTINUE 
50 CONTINUE 

CIRCRHO= (S*COS (PHI) *sumx) +(S*SIN(PHI) *sumy) + (C*sumz) 

RETURN 

END 


C*******FUNCTION CIRCPHI. 

C FUNCTION : CIRCPHI 

C DATE : 1 OCTOBER 1991 

C REVISED : 17 Feb 1992 

C PROGRAMMER : R.M. FRANCIS 

COMPLEX FUNCTION CIRCPHI(CNPHI,XP,AP,CNRHO,XR,AR,ARAD, 

* SCAN,PHI,RHB,ZHB) 

INTEGER CNPHI,CNRHO,M,N 

REAL XP(100),AP(100),XR(50),AR(50), 

* ARAD,SCAN,PHI,RHB,ZHB 
C LOCAL VARIABLES 

REAL Z,RZ,PHIPRIME,PI,S,C, 

* R,RP,X2,Y1,Y2,K,RR1,RR2 
COMPLEX J,JK,G1,G2,X1,C0N,CC, 

* sumx,sumy,sumz 
PI=3.141592654 

K=6.283185307 
RR1=(ARAD-0.)/2. 

RR2=(ARAD+0.)/2. 

J = (0.0,1.0) 

JK= (0.0,6.283185307) 
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CON=(0.0,-15.0) 

CC=C0N*RR1 

C OUTER LOOP: INTEGRATE OVER PHI. . . -PKPHKPI. 

C INNER LOOP: INTEGRATE OVER RHO... 0 <RHO<ANT 
RZ=RHB 
Z=ZHB 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 
surox=(0.0,0.0) 
sumy=(0.0,0.0) 
sumz=(0.0,0.0) 

DO 50 M=1,CNPHI 

PHIPRIME=PI*XP(M) 

DO 60 N=1,CNRHO 
RP=RR1*XR(N)+RR2 

R=SQRT(((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
Gl=(((K*R)**2)-l-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**5 
X1=JK*((RP*COS(PHIPRIME)*SIN(SCAN))-R) 

Yl=(RZ*COS(PHI))-(RP*COS(PHIPRIME)) 

X2=Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 

sumx=sumx+(CC*AP(M)*AR(N)*(G1+(X2*G2))*CEXP(X1)*RP) 
sumy=sumy+(CC*AP(M)*AR(N)*Y1*Y2*G2*CEXP(XI)*RP) 
sumz=suinz+(CC*AP(M) *AR(N) *Y1*Z*G2*CEXP(X1) *RP) 

60 CONTINUE 
50 CONTINUE 

CIRCPHI= (-SIN (PHI) *suxnx) 4- (COS (PHI) *sumy) 

RETURN 

END 


C*********SUBROUTINE TESTSPHERE 


C 

C 

C 

C 

C 


SUBROUTINE 

DATE 

PROGRAMMER 

COMMENTS 


TESTSPHERE 
21 FEBRUARY 1992 
R. M. FRANCIS 

GENERATES A SPHERE, WITH PLOTTED POINTS 
AT EQUAL INTERVALS IN THETA. 


SUBROUTINE TESTSPHERE(NP,ZH,RH,BASE,RS,ZP) 
INTEGER NP,I 

REAL ZH(400),RH(400),BASE,RS,ZP,SPRAD 

PI=3.1435926 

BK=2.*PI 

ZP=0.0 

WRITE(6,*)'ENTER AN EVEN NUMBER OF POINTS (NP ):' 
READ(5,*)NP 

WRITE(6,*)'ENTER SPHERE RADIUS' 

READ(5,*)SPRAD 
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BASE=SPRAD 
RS=SPRAD 
DO 1241 1=1,NP 

ANGLE=PI*FLOAT(1-1)/(2.*FLOAT(NP-1)) 

ZH(I)=BK*SPRAD*COS(ANGLE) 

RH(I)=BK*SPRAD*SIN(ANGLE) 

1241 CONTINUE 
RETURN 
END 

B. SAMPLE DATA FILE OPTLDBOR 

*** BOR RADIATION PATTERN FOR A CIRCULAR DISC RADIATOR USING 
GENEX *** 


BOR DIAMETER (WAVELENGTHS)=***** 
NUMBER OF GENERATING POINTS (NP)= 70 

SURFACE RADIUS50.00 2PRIME 0.00 


INDEX 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 
27 


NT NPHI 
4 60 


Z(I) 

RHO(I) 

50.000 

0.000 

49.987 

1.138 

49.948 

2.276 

49.883 

3.412 

49.793 

4.547 

49.676 

5.679 

49.534 

6.808 

49.366 

7.934 

49.173 

9.056 

48.954 

10.173 

48.710 

11.285 

48.440 

12.390 

48.146 

13.490 

47.826 

14.582 

47.482 

15.667 

47.113 

16.744 

46.720 

17.812 

46.302 

18.871 

45.861 

19.920 

45.395 

20.959 

44.906 

21.987 

44.394 

23.003 

43.859 

24.008 

43.301 

25.000 

42.721 

25.979 

42.119 

26.945 

41.494 

27.897 


SURF 

IMPED 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 

10000.00 

0.00 
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28 

40.848 

28.834 

10000.00 

0.00 

29 

40.182 

29.756 

10000.00 

0.00 

30 

39.494 

30.663 

10000.00 

0.00 

31 

38.786 

31.554 

10000.00 

0.00 

32 

38.057 

32.429 

10000.00 

0.00 

33 

37.309 

33.287 

10000.00 

0.00 

34 

36.542 

34.128 

10000.00 

0.00 

35 

35.755 

34.951 

10000.00 

0.00 

36 

34.951 

35.755 

10000.00 

0.00 

37 

34.128 

36.542 

10000.00 

0.00 

38 

33.287 

37.309 

10000.00 

0.00 

39 

32.429 

38.057 

10000.00 

0.00 

40 

31.554 

38.786 

10000.00 

0.00 

41 

30.663 

39.494 

10000.00 

0.00 

42 

29.756 

40.182 

10000.00 

0.00 

43 

28.834 

40.848 

10000.00 

0.00 

44 

27.897 

41.494 

10000.00 

0.00 

45 

26.945 

42.119 

10000.00 

0.00 

46 

25.979 

42.721 

10000.00 

0.00 

47 

25.000 

43.301 

10000.00 

0.00 

48 

24.008 

43.859 

10000.00 

0.00 

49 

23.003 

44.394 

10000.00 

0.00 

50 

21.987 

44.906 

10000.00 

0.00 

51 

20.959 

45.395 

10000.00 

0.00 

52 

19.920 

45.861 

10000.00 

0.00 

53 

18.871 

46.302 

10000.00 

0.00 

54 

17.812 

46.720 

10000.00 

0.00 

55 

16.744 

47.113 

10000.00 

0.00 

56 

15.667 

47.482 

10000.00 

0.00 

57 

14.582 

47.826 

10000.00 

0.00 

58 

13.490 

48.146 

10000.00 

0.00 

59 

12.390 

48.440 

10000.00 

0.00 

60 

11.285 

48.710 

10000.00 

0.00 

61 

10.173 

48.954 

10000.00 

0.00 

62 

9.056 

49.173 

10000.00 

0.00 

63 

7.934 

49.366 

10000.00 

0.00 

64 

6.808 

49.534 

10000.00 

0.00 

65 

5.679 

49.676 

10000.00 

0.00 

66 

4.547 

49.793 

10000.00 

0.00 

67 

3.412 

49.883 

10000.00 

0.00 

68 

2.276 

49.948 

10000.00 

0.00 

69 

1.138 

49.987 

10000.00 

0.00 

70 

0.000 

50.000 

10000.00 

0.00 


PHI OF RECEIVER (DEG)= 0.00 

MAXIMUM FIELD VALUE (V/M)= 0.16659E+03 
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ANGLE 


co- 

■POL 



X-POL 

(DEG) 

(VOLTS) 

(DEG) 

(DB-REL) 

(VOLTS) 

(DEG) 

(DB-REL) 

0.00 

6.9271 

-104.51 

-27.62 

0.0000 

-92.95 

-100.00 

1.30 

3.6864 

-106.76 

-33.10 

1.1253 

131.93 

-43.41 

2.61 

2.6516 

-50.40 

-35.96 

0.7957 

163.87 

-46.42 

3.91 

4.8350 

-14.07 

-30.74 

0.5577 

-172.05 

-49.50 

5.22 

6.0401 

-5.71 

-28.81 

0.2076 

-139.66 

-58.09 

6.52 

4.7465 

-2.93 

-30.91 

0.2741 

5.44 

-55.68 

7.83 

1.3824 

-4.09 

-41.62 

0.5842 

30.18 

-49.10 

9.13 

2.8075 

-174.09 

-35.47 

0.7556 

57.97 

-46.87 

10.43 

5.8878 

-174.09 

-29.03 

0.9056 

82.84 

-45.29 

11.74 

6.3663 

-172.38 

-28.36 

0.9042 

101.87 

-45.31 

13.04 

3.9668 

-170.31 

-32.46 

0.6473 

118.37 

-48.21 

14.35 

0.9093 

7.83 

-45.26 

0.2436 

-130.52 

-56.70 

15.65 

6.0107 

6.17 

-28.85 

1.0247 

-96.02 

-44.22 

16.96 

8.2390 

4.51 

-26.12 

1.2154 

-95.19 

-42.74 

18.26 

6.7229 

7.30 

-27.88 

0.9426 

-80.30 

-44.95 

19.57 

2.3749 

15.98 

-36.92 

0.4040 

-54.39 

-52.31 

20.87 

3.7384 

-170.33 

-32.98 

0.5796 

72.59 

-49.17 

22.17 

9.6622 

-170.82 

-24.73 

1.3618 

88.39 

-41.75 

23.48 

10.7713 

-175.50 

-23.79 

1.5721 

81.57 

-40.50 

24.78 

6.6302 

-160.18 

-28.00 

0.6277 

115.25 

-48.48 

26.09 

3.5922 

-149.07 

-33.33 

0.8164 

137.53 

-46.20 

27.39 

7.2792 

27.16 

-27.19 

0.9508 

-75.15 

-44.87 

28.70 

13.9956 

3.03 

-21.51 

2.1212 

-107.94 

-37.90 

30.00 

12.6783 

26.06 

-22.37 

1.1997 

-71.27 

-42.85 

31.30 

14.4138 

12.62 

-21.26 

2.0480 

-86.02 

-38.21 

32.61 

1.4466 

93.02 

-41.23 

0.5246 

44.99 

-50.04 

33.91 

5.0101 

-168.15 

-30.44 

0.7135 

42.79 

-47.37 

35.22 

18.2410 

-161.85 

-19.21 

2.3460 

82.73 

-37.03 

36.52 

22.8067 

-158.28 

-17.27 

2.5909 

93.03 

-36.16 

37.83 

17.8045 

-160.35 

-19.42 

2.0096 

80.80 

-38.37 

39.13 

18.0824 

-146.48 

-19.29 

2.4138 

124.03 

-36.78 

40.43 

9.1689 

3.14 

-25.19 

1.7552 

-103.63 

-39.55 

41.74 

13.8129 

38.91 

-21.63 

0.9059 

-103.47 

-45.29 

43.04 

30.3979 

9.45 

-14.78 

4.3773 

-120.67 

-31.61 

44.35 

47.4641 

29.97 

-10.91 

5.2557 

-84.81 

-30.02 

45.65 

36.1417 

33.25 

-13.27 

3.3511 

-87.55 

-33.93 

46.96 

31.6252 

32.13 

-14.43 

3.2704 

-81.37 

-34.14 

48.26 

14.5891 

62.48 

-21.15 

2.0676 

-25.23 

-38.12 

49.57 

18.9270 

-170.31 

-18.89 

2.8516 

41.68 

-35.33 

50.87 

45.8795 

-157.68 

-11.20 

5.2496 

54.10 

-30.03 

52.17 

78.0162 

-154.23 

-6.59 

8.1824 

59.69 

-26.18 

53.48 

106.099 

-152.60 

-3.92 

10.9266 

61.19 

-23.66 

54.78 

136.447 

-152.71 

-1.73 

13.7615 

65.23 

-21.66 

56.09 

160.822 

-147.27 

-0.31 

14.3460 

71.03 

-21.30 

57.39 

161.358 

-145.24 

-0.28 

13.9190 

68.77 

-21.56 

58.70 

166.587 

-146.15 

0.00 

14.2015 

67.31 

-21.39 

60.00 

162.565 

-144.37 

-0.21 

13.4435 

67.27 

-21.86 

61.30 

153.157 

-144.83 

-0.73 

12.5463 

66.68 

-22.46 

62.61 

143.692 

-144.21 

-1.28 

11.3471 

67.74 

-23.34 
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65.22 

108.369 

-140.84 

-3.73 

7.6142 

65.73 

-26.80 

66.52 

85.0397 

-139.06 

-5.84 

5.6224 

60.07 

-29.43 

67.83 

67.7272 

-140.66 

-7.82 

4.2600 

58.69 

-31.84 

69.13 

59.7937 

-143.81 

-8.90 

3.8263 

67.63 

-32.78 

70.43 

47.3499 

-145.20 

-10.93 

3.1370 

62.89 

-34.50 

71.74 

30.9244 

-145.07 

-14.63 

2.1385 

43.54 

-37.83 

73.04 

19.4698 

-143.32 

-18.65 

1.2247 

37.25 

-42.67 

74.35 

12.7001 

-149.46 

-22.36 

0.8089 

32.28 

-46.28 

75.65 

5.2649 

-153.77 

-30.01 

0.8572 

-15.23 

-45.77 

76.96 

8.7952 

-119.58 

-25.55 

0.7979 

-159.71 

-46.39 

78.26 

12.4865 

173.28 

-22.50 

1.4228 

77.84 

-41.37 

79.57 

11.0782 

91.40 

-23.54 

1.4678 

-11.88 

-41.10 

80.87 

3.7571 

-1.39 

-32 . < 

1.0531 

-129.37 

-43.98 

82.17 

6.6444 

103.97 

-27.9c 

0604 

30.52 

-43.92 

83.48 

6.3822 

-4.11 

-28.33 

68 

-128.96 

-43.87 

84.78 

3.0249 

-158.22 

-34.82 

3 4 

104.94 

-54.94 

86.09 

1.8402 

75.64 

-39.14 

0.2566 

-5.90 

-56.25 

87.39 

3.9405 

-150.16 

-32.52 

0.8277 

166.07 

-46.08 

88.70 

7.3638 

110.06 

-27.09 

1.1926 

21.33 

-42.90 

90.00 

7.6190 

-5.04 

-26.79 

1.2270 

-108.00 

-42.66 


C. SOURCE CODE FOR GENPLT.M 

%PROGRAM : genplt.m 

%DATE : 15 February 1992 

%PROGRAMMER: R.M. FRANCIS 

%Comments : This prc-ram is ? plotting routine for the data 
% points generated : -•••- ? - -rv rrar 7..BGRMM.F and the 

% program TEST.M. tne runs of test.m and 

% LDBORMM.F. 


load bang.m 
load bxpole.m 
load bcpole.m 
load ctest 
for i=l:np 

if bcpole(i)<=0.0 
bcpole(i)=eps; 

end 

if bxpole(i)<=0.0 
bxpole(i)=eps; 

end 

el(i)=10*log(bcpole(i)) ; 
ep(i)=10*log(bxpole(i)); 
zl(i)=10*log(z(i)) ; 
if zl(i)<-50 
zl(i)=-50; 

end 
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if el(i)<-50 
el(i)=-50; 

end 

if ep(i)<-50 
ep(i)=-50; 

end 

end 

subplot(221),plot(deg,z,bang,bcoole,'.') 
title('TEST OF LDBORMM. F') ,xlabel(' degrees') 
ylabel('magnitude') 

subplot (222) , plot (deg, zl,b_ng,el, '.') 
title('LOG PLOT (db)'),xlabel('degrees') 
ylabel('magnitude') 

text(.1,.4,'analytical (solid), LDBORMM.F(.)','sc') 
text(.1,.3,'PARAMETERS:','sc') 
text(.1,.25,'number of points = 90','sc') 
text(.5,.25,'scan angle = 60 deg.','sc') 
text(.5,.15,'CNRHO = 20','sc') 
text(.5,.2,'CNPHI = 60','sc') 
text(.1,.2,'NPHI = 100','sc') 
text(.l,.15,'NT = 10','sc') 
text(.1,.1,'distance = 50','sc') 
text(.5,.1,'aperture radius = 5','sc') 
text(.1,.05,'complex impedance = (10000.0, 0.0)','sc') 
pause 
%print 
clg 

load etscat.m 
load epscat.m 
load etf.m 
load epf.m 

subplot(221),plot(bang,etscat) 
title('ETSCAT') 

subplot(222),plot(bang,epscat) 
title('EPSCAT') 
subplot(223),plot(bang,etf) 
title('ETF') 

subplot(224),plot(bang,epf) 
title('EPF') 

%print 
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APPENDIX B 

TABLE B.l. VARIABLES IN SUBROUTINE GENEXX 


VARIABLE 

TYPE 

DESCRIPTION 

RANGE 

NM 

integer 

Number of modes in 0 

1< NM < oo 

NP 

integer 

Even number of points on 
surface of the BOR 

2< NP <400 

NT 

integer 

Even number of points of 
integration in T on BOR 

2< NT <4 

NPHI 

integer 

Even number of points of 
integration in 0 on BOR 

2<NPHI<100 

CNRHO 

integer 

Even number of points of 
integration ir. p 1 (source) 

2<CNRHO<50 

CNPHI 

integer 

Even number of points of 
integration in 0* (source) 

2<CNPHI<100 

XT, X, XR, 
XP 

real 

Abscissa of the function (x n ) 
in T,0,p ' and 0' 

-i<x n <i 

AT,A,AR, 
AP 

real 

Weights of the series 
terms(for x Q ) approximating 
the integral 

0<o) n <l 

SCAN 

real 

angular displacement (in 
degrees) of source boresight 
and Z axis 

-9O<0 s <9O 

PHI 

real 

orientation in <p of receiver 

O<0<18O 

ARAD 

real 

dimensionless radius of 
source antenna (a/A) 

0<p ' < oo 

RH(i) 

real 

the product of the wave 
number (k) and the radius of 
the BOR, units in radians 

0< r < oo 

3< i < 399 

ZH(i) 

real 

the product of the wave 
number (k) and the distance z 
on the BOR axis, units in 
radians 

0< z < oo 

3< i < 399 

R(i) 

complex 

the excitation vector 
i =2NP-3 

max 

- oo< R(i) <oo 

4< i < 1600 
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APPENDIX C 

A. PROGRAM TO GENERATE DATA FILES FOR GAUSS QUADRATURE 
INTEGRATION 

C PROGRAM : GAUS.F 

C DATE : 12 SEPTEMBER 1991 

C PROGRAMMER : R. FRANCIS, D. JENN 

C PROGRAM TO COMPUTE WEIGHTS AND ABSCISSAS FOR GAUSSIAN INT. 
CHARACTER*? FNAME 
DIMENSION X(500),A(500) 

WRITE(6,*) 'ENTER FILE NAME (gaus-)' 

READ(5,*)FNAME 

WRITE(6,*) 'ENTER AN EVEN NUMBER (max pts is 500)' 

READ(5,*) N 

OPEN(2,FILE=FNAME) 

WRITE(2,*) N 

CALL GAUSLEG(-1.,1.,X,A,N) 

DO 10 1=1,N 

10 WRITE(2,*) X(I),A(I) 

STOP 

END 

SUBROUTINE GAUSLEG(XI,X2,X,W,N) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*4 X1,X2,X(N),W(N) 

PARAMETER (EPS=3.D-14) 

M=(N+l)/2. 

XM=.5D0*(X2+X1) 

XL=.5D0*(X2-X1) 

DO 12 1=1,M 

Z=COS(3•1415926535D0*(I -.25D0)/(N+.5D0)) 

1 CONTINUE 

Pl=l.DO 
P2=0.DO 
DO 11 J=1,N 
P3=P2 
P2=P1 

Pl=((2.D0*J-1.D0)*Z*P2-(J-l.DO)*P3)/J 

11 CONTINUE 

PP=N*(Z*P1-P2)/(Z*Z-1.DO) 

Z1=Z 

Z=Z1-P1/PP 

IF(ABS(Z-Z1).GT.EPS) GO TO 1 
X(I)=XM-XL*Z 
X(N+l-I)=XM+XL*Z 
W(I)=2.D0*XL/((l.D0-Z*Z)*PP*PP) 

W(N+l-I)=W(I) 

12 CONTINUE 
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RETURN 

END 


B. SAMPLE DATA FILE GENERATED BY GAPS.F 

The following is an example of the external sequential 
data file generated by the program GAUS.F. The first row is an 
even integer (N) . The data in rows two through N are the 
abscissa and weights for the series approximation of the 
integral calculation. 

Sample data file: qaus20 _ 


20 

-0.993129 

1.76140E-02 

-0.963972 

4.06014E-02 

-0.912234 

6.26720E-02 

-0.839117 

8.32767E-02 

-0.746332 

1.01930E-01 

-0.636054 

0.118195 

-0.510867 

0.131689 

-0.373706 

0.142096 

-0.227786 

0.149173 

-7.65265E-02 

0.152753 

7.65265E-02 

0.152753 

0.227786 

0.149173 

0.373706 

0.142096 

0.510867 

0.131689 

0.636054 

0.118195 

0.746332 

1.01930E-01 

0.839117 

8.32767E-02 

0.912234 

6.26720E-02 

0.963972 

4.06014E-02 

0.993129 

1.76140E-02 
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APPENDIX D 


A. ARGUMENTS: CIRCTHETA. CIRCRHO AND CIRCPHI 
TABLE D.l. VARIABLE LIST FOR FUNCTION SUBPROGRAMS. 


VARIABLE 

TYPE 

DESCRIPTION 

RANGE 

CNPHI 

integer 

see Table B.l. 

same 

CNRHO 

integer 

see Table B.l 

same 

XP, XR 

real 

see Table B.l 

same 

AP, AR 

real 

see Table B.1 

same 

ARAD 

real 

see Table B.l 

same 

SCAN 

real 

see Table B.l 

same 

PHI 

real 

see Table B.l 

same 

RHB 

real 

same as RH(i) of 

Table B.l 

same 

ZHB 

real 

same as ZH(i) of 

Table B.l 

same 

circtheta 

circrho 

circphi 

complex 

E e , E p and Ep 
calculated at the 
point on the BOR 
surface defined by 

RHB,ZHB and PHI 

all 

complex 

numbers 


B. TEST PROGRAM: CIRCSUB.F 

The following program was used to test for convergence the 
numerical methods used to calculate the electric field 
components for a uniformly illuminated circular aperture. 
Before executing the main program LDBORMM.F, it is desirable 
to determine the minimum number of terms required (CNPHI and 
CNRHO) for the series approximation of the source field to 
converge to the analytical solution of Equation (3.26). The 
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following programs are provided to test for the number of 
terms required for convergence. 


The programs CIRCSUB.F, TEST.M and SCANPLT.M are executed 
in the same directory. To minimize the number of files in the 
main directory, maintain these programs in a subdirectory 
named "circtest”. CIRCSUB.F creates the external sequential 
data files ang.m, etheta.m, erho.m and ephi.m. These data 
files are vectors of length NP, containing the angle (0) and 
the magnitude of E 0 , E p and E^. The matlab program TEST.M 
calculates the analytic solution for E fl in accordance Equation 
(3.26) and writes the data to a file CTEST.MAT. The matlab 
program SCANPLT.M loads the data files and plots the results. 

1. Source Code for CIRCSUB.F 


C TESTPROGRAM : CIRCSUB.F 

C DATE : 1 OCTOBER 1991 

C REVISED : 12 JANUARY 1992 

C PROGRAMMER : R.M. FRANCIS 

C COMMENTS : THIS SUBROUTINE CALCULATES THE ELECTRIC FIELD 
C FOR A CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE 
C BOUNDARY DEFINED BY AN OGIVE. THE APERTURE IS LOCATED AT Z 
C = 0, AND IS BORE SIGHTED ON THE Z AXIS. THE SUBROUTINE 
C OGIVE IS THE SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC 
C TO PERFORM COMPUTATIONS. 

C ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 


C********* PROGRAM CIRC 

C DECLARE VARIABLES INTEGER NP,NPHI,NRHO,I,M,N 
C THE CHARACTER STRING PHIPTS/RHOPTS WILL ACCESS THE 
C APPROPIATE FILE WITH THE X'S AND COEEFFICIENTS, TO PERFORM 
C INTEGRATION BY GAUSSIAN QUADRATURE. EG. GAUS10 YEILDS, 10 
C POINTS;GAUS20, 20PTS, ETC. 

CHARACTER*15 PHIPTS,RHOPTS,GP,GR 

REAL XP(400),AP(400),XR(100),AR(100),S,C,P,PI,SPRAD, 

* SCAN,Z,RZ,PHI,PHIPRIME,DEG,SA,K,ET,ER,EP, 
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* R, RP,X2,Y1,Y2,ARAD,RH(400),ZH(400),angle(400) 
COMPLEX J,JK,G1,G2,X1,C0N,CC, 

* sumx, suiny, sumz, 

* ETHETA(400),ERHO(400),EPHI(400) 

PI=3.141592654 

J = (0.0,1.0) 

JK= (0.0,6.283185307) 

K=6.283185307 
CON=(0.0,-15.0) 

WRITE(6,*)'ENTER THE RADIUS OF THE ANTENNA,(in wavelengths)' 
READ(5,*) ARAD 
ANT=ARAD/2. 

WRITE(6,*)•ENTER PHI,(OBSERVATION) IN DEGREES' 

READ(5,*)P 
PHI=P*PI/180. 

WRITE(6,*)'ENTER SCAN ANGLE IN DEGREES' 

READ(5,*)SA 
SCAN=-SA*PI/180. 

WRITE(6,*)'ENTER PHIPTS STRING: "gaus### ",(max 400)' 
READ(5,*)GP 
PHIPTS='gaus/'//GP 

WRITE(6,*)'ENTER RHOPTS STRING: "gaus### ",(max 100)' 
READ(5,*)GR 
RHOPTS='gaus/•//GR 
CALL OGIVE(NP,ZH,RH) 
make loop to enter test sphere. 

write(6,*)'enter the sphere radius' 
read(5,*)SPRAD 

write(6,*)'enter the number of points' 
read(5,*)NP 
do 888 1=1,NP 

angle(I)=PI*(1-1)/(2*(NP-1)) 

RH(I)=SPRAD*SIN(angle(I))*K 
ZH(I)=SPRAD*COS(angle(I))*K 
write(6,*)angle(I),RH(I),ZH(I) 

888 continue 

OPEN(2,FILE=PHIPTS) 

OPEN(3,FILE=RHOPTS) 

OPEN(4,FILE='etheta.m') 

OPEN(8,FILE='erho.m') 

OPEN(9,FILE='ephi.m') 

OPEN(10,FILE='ang.m') 

READ(2,*)NPHI 

READ(3,*) NRHO 

DO 20 M=1,NPHI 

READ(2,*,END=2 0)XP(M),AP(M) 

20 CONTINUE 

DO 30 N=1,NRHO 

READ(3,*,END=30)XR(N),AR(N) 

30 CONTINUE 

C OUTER LOOP : SOLVES FOR THE SPHERICAL COMPONENTS OF E AT 
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C EACH POINT ON THE SURFACE OF THE OGIVE. 

C MIDDLE LOOP: INTEGRATE OVER PHI...-PI<PHI<PI. 

C INNER LOOP: INTEGRATE OVER RHO... 0 <RHO<ANT 
DO 40 1=1,NP 

Z=ZH(I)/K 

RZ=RH(I)/K 

S=RZ/SQRT(RZ**2+Z**2) 

C= Z/SQRT(RZ**2+Z**2) 

DEG=angle(I)*180/PI 
sumx=(0.0,0.0) 
sumy=(0.0,0.0) 
sumz=(0.0,0.0) 

DO 50 M=1,NPHI 

PHIPRIME=PI*XP(M) 

DO 60 N=1,NRHO 
RP=(ANT*XR(N))+ANT 

R=SQRT((RZ-RP)**2+(Z**2)+(4*RZ*RP*SIN((PHI-PHIPRIME)/2)**2)) 
Gl=(((K*R)**2)-l-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**5 
X1=JK*((RP*COS(PHIPRIME)*SIN(SCAN))-R) 

Yl=(RZ*COS(PHI))-(RP*COS(PHIPRIME)) 

X2=Y1**2 

Y2=(RZ*SIN(PHI))-(RP*SIN(PHIPRIME)) 

CC=CON*ANT 

sumx=surax+(CC*AP(M)*AR(N)*(G1+(X2 *G2))*CEXP(XI)*RP) 
sumy=sumy+(CC*AP(M)*AR(N)*Y1*Y2*G2*CEXP(XI)*RP) 
sumz-sumz+(CC*AP(M)*AR(N)*Y1*Z*G2*CEXP(X1)*RP) 

60 CONTINUE 
50 CONTINUE 

C CALCULATE THE SPHERICAL COMPONENTS AT THE SURFACE. 

ETHETA(I)=((C*COS(PHI)*sumx)+(C*SIN(PHI)*sumy)-(S*sumz)) 

ERHO(I) =((S*COS(PHI)*sumx)+(S*SIN(PHI)*sumy)+(C*sumz)) 
EPHI(I) =((-SIN(PHI)*sumx)+(COS(PHI)*sumy)) 

ET=CABS(ETHETA(I)) 

ER=CABS(ERHO(I)) 

EP=CABS(EPHI(I)) 

WRITE(4,*)ET 
WRITE(8,*)ER 
WRITE(9,*)EP 
WRITE(10,*)DEG 
40 CONTINUE 
STOP 
END 


4 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OGIVE 

4 SEPTEMBER 1991 
R.M. FRANCIS 
13 JANUARY 1992 


SUBROUTINE 
DATE 

PROGRAMMER 
REVISED 
COMMENTS 

THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION 
(BOR) IN THE FORM OF AN OGIVE. 

DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 

ZH = Z CO-ORDINATE * 2*PI 
RH = RADIUS *2*PI 


C***************** 

SUBROUTINE OGIVE(NP,ZH,RH) 

C NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 
C ZP = ZPRIME, THE POSITION ON Z WHERE THE RADIUS OF CURVATURE 
C STARTS 

C B = BASE RADIUS 

C RS = RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 
INTEGER I,NP 

REAL RH(400),ZH(400),ZP,BASE,RS,ZCOORD,RADIUS,PI 
PI= 3.141592654 

C INPUT THE VARIABLES FOR THE OGIVE,ZP,B,RS,NP 

WRITE(6,*)'ENTER AN EVEN NUMBER OF POINTS, MAX IS 400' 
READ(5,*)NP 

WRITE(6,*) 'ENTER SURFACE CURVATURE (wavelengths)' 

READ(5,*) RS 

WRITE(6,*) 'ENTER ZPRIME,WHERE CURVATURE STARTS (wavelengths)' 
READ(5,*) ZP 

WRITE(6,*) 'ENTER BASE RADIUS (wavelengths)' 

READ(5,*) BASE 
C PERFORM CALCULATIONS 

ZMAX= SQRT((2*BASE*RS)-BASE**2) + ZP 
DZ= ZMAX/(float(NP)-l.) 

DO 10 1=1,NP 

ZCOORD= (float(I)-1.)*DZ 
ZH(I)=2.*PI*ZCOORD 

IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 

ELSE 

RADIUS=SQRT(RS**2-(ZCOORD-ZP)**2)+(BASE-RS) 
ENDIF 

RH(I)=2.*PI*RADIUS 
10 CONTINUE 
RETURN 
END 
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2. Source Code for TEST.M 


% PROGRAM : TEST.M 

% DATE : 3 OCTOBER 1992 

% REVISED : 15 FEBRUARY 1992 

%Program TEST.M— generate analytic soln. for circular 
% aperture. 

%saves the data to "mat file", ctest. 
npts=input('enter the number of points:') 

R=input('enter the distance (R):') 

r=input('enter the radius (in wavelengths): ') 

scan=input('enter the scan angle(degrees)') 

sa=scan*pi/180; 

eta= 120*pi; 

k=2*pi; 

j=sqrt(-1); 

nmax=2*(npts-l) ; 

cl=j *k*eta*r A 2*exp(-j *k*R)/R; 

% calculate the solution from 0 to 90 degrees.******* 
for i=0:npts-l 
ang(i+l)=pi*i/nmax; 
deg(i+l)=90*i/(npts-l); 
fac=0.5*cos(ang(i+l)); 
if sin(ang(i+l))-sin(sa)== 0.0 
a(i+l)=eps; 

else 

a(i+l)=2*pi*r*(sin(ang(i+l))-sin(sa)); 

end 


end 


z(i+l)=fac*abs(cl*bessel(1,a(i+l))/a(i+l)); 


np=npts; 

save ctest deg z np 


3. Data file 

0. OOOOOOOe-^OO 
4.7368421e+00 
9.4736842e+00 
1.4210526e+01 
1.8947368e+01 
2.3684211e+01 
2.8421053e+01 
3.3157895e+01 
3.78947 37e+01 
4.2631579e+01 
4.7368421e+01 
5.2105263e+01 
5.6842105e+01 
6.1578947e+01 
6.6315789e+01 
7.1052632e+01 
7.5789474e+01 


CTEST.MAT 

5.9217626e+01 
5.7051528e+01 
5.0941450e+01 
4.1944793e+01 
3.1506384e+01 
2.1103366e+01 
1.1933350e+01 
4.7301405e+00 
2.6730656e-01 
3.2224083e+00 
4.5450714e+00 
4.7331488e+00 
4.2520406e+00 
3.4672834e+00 
2.6244918e+00 
1.8607876e+00 
1.2307921e+00 


20 
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8.0526316e+01 7.3451071e-01 
8.5263158e+01 3.4018017e-01 
9.0000000e+01 2.4512433e-16 

4. Source code for SCANPLT.M 

%PROGRAM : scanplt.m 

%DATE : 15 February 1992 

%Comments : This program is a plotting routine for the data 
% points generated by the test program CIRCSUB.F and the 
% program TEST.M. Enter parameters iaw the runs of 
% TEST.M and CIRCSUB.F. 
load ang.m 
load etheta.m 
load ctest 

subplot(211),plot(deg,z,ang,etheta,'.') 

title('TEST OF SUBROUTINE CIRCTHETA'),xlabel('degrees') 
ylabel( 1 magnitude') 

text(.1,.4analytical (solid), CIRCTHETA (.)','sc') 

text(.1,.3, 1 PARAMETERS:','sc') 

text(.1,.25,'number of points = 90','sc') 

text(.5,.25,'scan angle = 60 deg.','sc') 

text(.1,.2, 'CNRHO * 20','SC' ) 

text(.5,.2,'CNPHI = 60','sc') 

text(.1,.15,'distance = 2500','sc') 

text(.5,.15,'aperture radius = 5','sc') 
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APPENDIX E. MATH NOTES 
A. PLOTS GENERATED BY SUBROUTINE OGIVE 

Plots generated by the subroutine OGIVE do not have equal 
spacing in the angle 0. As the angle increases the points are 
spaced more closely. Equation 2.1 defines r(z). It follows 
that 

0 = tan' 1 | 1 j (E. l) 

as z approaches zero, 0 approaches 90 degrees. The rate at 
which 0 approaches 90 degrees decreases as z decreases since 


d 0 . z 2 
d z z 2 * z (z) 2 

clearly demonstrates in the limit 


limit 
z- 0 


d 0 
d z 


= 0 


(E. 2) 


(E.3) 


Where the subroutine OGIVE generates the points RH and ZH, the 
density of plotted values increases with angle 0. 
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B. PLOTS GENERATED BY SUBROUTINE TESTSPHERE 

In the subroutine TESTSPHERE, the variables RH and ZH are 
defined in terms of a sphere radius (SPRAD) and the angle 0 . 
The following pseudocode applies 

e (i) = e scart ♦ a 0 

RH(i) = SPRAD * sin (0 ( i ) ) ( E - 4 > 

ZHU) = SPRAD * cos ( 0 ( i ) ) 

where A0 is a constant given by 

and NP is the number of points on the BOR. 


A 0 = _ Q ct ? £t (E.5) 

NP - 1 

C. NEAR FIELD EQUATIONS FOR A CIRCULAR APERTURE 

The Equations 3.13 to 3.22 were adapted from Equations ( 6 - 
108a) to (6-108c) on page 283 of reference 5. To implement the 
near field integrals for a circular aperture bounded by a BOR 
the cartesian representation of the equations in reference ?? 
were converted to polar coordinates as follows 

x' = p' cos <J>' 

y' = p' sin <t>' (E. 6) 

ds 1 = dx' dy' = p 1 dp' d4> ; 
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Where the primed notation denotes a source coordinate. 
Straight forward substitution of Equation (E.6) into the 
equations of reference 5 yeild Equations 3.13 to 3.22. 

The cartesian components of the electric field E(x,y,z) 
are transformed to spherical components by 

E e [cos (0) cos (<j>) cos (6) sin (<$>) -sin (0)1 E x 

E p = sin (0) cos (<f>) sin(0) sin(<{>) cos(0) E y ( E * 7 > 

Em. -sin(4>) cos($) 0 e, 

▼J 2. 



Figure E.l. Polar Geometry for Conversion of E(x,y,z) to 
Spherical Components,E(R,0, <p) . 

Figure E.l. illustrates the geometry. 
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